Methods and computer software for detecting splice variants

ABSTRACT

Methods and software products for analysis of alternative splicing are disclosed. In general the methods involve normalizing probe set or exon intensity to an expression level measurement of the gene. The methods may be used to identify tissue-specific alternative splicing events.

RELATED APPLICATIONS

The present application claims priority to U.S. Provisional ApplicationNo. 60/722,742, filed Sep. 30, 2005, the disclosure of which isincorporated herein by reference in its entirety.

FIELD OF THE INVENTION

The invention is related to computer methods to identify differentialsplicing events. Specifically, this invention provides methods, computersoftware products and systems for analysis of microarray data toidentify differentially spliced exons.

BACKGROUND OF THE INVENTION

Recent genome-wide analysis of alternative splicing indicates that alarge portion of human genes, probably more than half, have alternativesplice forms. Alternative splicing provides the cell with a mechanism togenerate multiple gene products from the same transcript, adding to thefunctional complexity of the genome. Regulated alternative splicing maybe used to create different proteins under different circumstances,allowing production of functionally related but distinct proteins andthus expanding the protein-coding potential of genes and genomes.

The identities of the genes that are being expressed in a biologicalsample at any given time and the amount of expression of those genesprovide a gene expression profile for that sample. The gene expressionprofile is an indication of the status of that sample. For example,different tissue types will have different gene expression profilesreflecting the expression of different genes and differences in thespliced forms of individual genes. Differences in expression profile mayalso be observed between samples from the same tissue type when onesample is diseased. High-throughput methods to analyze and detectexpression of alternative splice forms, characterization of alternativesplicing, and regulation of alternative splicing are an importantresearch focus.

SUMMARY OF THE INVENTION

Computer software products and methods for analysis of alternativesplicing are disclosed.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1: Schematic of Exon Array data analysis workflow.

FIG. 2: ROC curve for Colon samples (Normal vs. Tumor) using MIDAS

FIG. 3: Histogram of exons per Constitutive Core Gene

FIG. 4: Comparison of Empirical p-Values with Theoretical p-Values

FIG. 5: The log-scale error in approximating background by threetechniques: global, gc-bin averages using BGP, and individual mismatches(MM) is shown here for 8,103 GC-bin probes.

FIG. 6: Despite the similarity in the boxplots, mismatched probes (MM)are still a better approximation in mean-squared error (log-scale).However, since there is some loss of signal due to MM hybridization tothe desired target, this extra variation does not greatly impact signalaccuracy.

FIG. 7: (Rule 2e) Extended annotations that do not directly overlap anycore annotations are not added to those associated clusters. Insteadthey are merged to form a cluster of their own.

FIG. 8: (Rule 2a) The arrows indicate a shared splice site between thetwo extended annotations. The bottommost extended annotation is joinedwith the cluster that it shares a splice site with.

FIG. 9: (Rule 2c) The extended annotation is broken up into underlyingexons because it overlaps two gene clusters from a higher confidencelevel.

FIG. 10: Transcript annotations from different confidence levels aremerged to form a gene annotation. The regions of the gene annotation canbe labeled according to the highest level confidence transcript thatsupports that region.

FIG. 11: Probesets are labeled with a confidence ranking according tothe confidence level region of the overlapping gene. Probesets that fallwithin multiple genes are labeled ‘ambiguous’, unless the probeset fallswithin exactly one core region of a gene; then it is labeled ‘core’.Probesets that overlap confidence region boundaries are labeled with thelower confidence level. Probesets that do not fall within any genes arelabeled ‘free’.

FIG. 12: Probesets mapping to non-core single exon annotations that fallwithin the intron of one other gene are grouped with that gene and giventhe label ‘bounded’.

FIG. 13: An annotation is evidence for a probeset if the probeset liescompletely within one of the exons of the annotation.

FIG. 14: Correlation with original estimates as low-signal decoy probesets were added. As more low-signal GENSCAN Suboptimal decoy probe setswere added to the gold standard probe sets, the correlation with theoriginal estimates decreases. However, the robustness of PLIER isevident as even when adding 10 unrelated probesets the correlationremains high.

FIG. 15: Correlation with original Plier estimates as high-signal decoyprobe sets were added. Including decoy probe sets from full length mRNAsdecreases the correlation with the original estimates. PLIER is not asresistant to higher-signal decoy probe sets as it is to low-signal decoyprobe sets (FIG. 16).

FIG. 16: Signal ROC Plot for Differentiating Between Spike-inConcentrations of 0.71 pM and 0.36 pM. With no decoy probe sets, theGeneChip Exon Array differentiates well between spike-ins at differentconcentrations and those added at the same concentration.

FIG. 17: Signal ROC Plot of Spike-in Concentration (0.71 pM vs 0.36 pM)Differentiation with 20 Decoy Probe Sets. The addition of low-signalprobe sets decreases the ability of the PLIER expression estimates todifferentiation between spike-in concentration pools.

FIG. 18: Correlation with original set as low-intensity decoys are addedusing an IterPLIER. The IterPLIER expression estimates are bettercorrelated to the gold standard PLIER expression estimates than the rawPLIER expression estimates.

FIG. 19: Correlation with original set as high-intensity decoys areadded using IterPLIER. The IterPLIER expression estimates are bettercorrelated to the gold standard PLIER expression estimates than the rawPLIER expression estimates.

FIG. 20: IterPLIER Improves Signal ROC performance in Presence of 20Decoy Probe Sets (0.71 pM vs 0.36 pM). The ability of IterPLIERexpression estimates to differentiate between spike-in concentrations isimproved compared to raw PLIER expression estimates.

FIG. 21: Boxplots of log 2 cell intensities for a set of 20 arrays.

FIG. 22: Histograms of log 2 cell intensities for the first 8 arrays inthe set of 20.

FIG. 23: Boxplots of M values for a set of arrays: M=relative log 2intensity=relative log 2 intensity=log 2 intensity for array−median log2 intensity

FIG. 24: Boxplots of M values for a set of arrays: M=relative log 2plier signal

FIG. 25: Bar chart of mad(plier residuals)

DETAILED DESCRIPTION OF THE INVENTION a) General

The present invention has many preferred embodiments and relies on manypatents, applications and other references for details known to those ofthe art. Therefore, when a patent, application, or other reference iscited or repeated below, it should be understood that it is incorporatedby reference in its entirety for all purposes as well as for theproposition that is recited.

As used in this application, the singular form “a,” “an,” and “the”include plural references unless the context clearly dictates otherwise.For example, the term “an agent” includes a plurality of agents,including mixtures thereof.

An individual is not limited to a human being, but may also be otherorganisms including, but not limited to, mammals, plants, bacteria, orcells derived from any of the above.

Throughout this disclosure, various aspects of this invention can bepresented in a range format. It should be understood that thedescription in range format is merely for convenience and brevity andshould not be construed as an inflexible limitation on the scope of theinvention. Accordingly, the description of a range should be consideredto have specifically disclosed all the possible subranges as well asindividual numerical values within that range. For example, descriptionof a range such as from 1 to 6 should be considered to have specificallydisclosed subranges such as from 1 to 3, from 1 to 4, from 1 to 5, from2 to 4, from 2 to 6, from 3 to 6 etc., as well as individual numberswithin that range, for example, 1, 2, 3, 4, 5, and 6. This appliesregardless of the breadth of the range.

The practice of the present invention may employ, unless otherwiseindicated, conventional techniques and descriptions of organicchemistry, polymer technology, molecular biology (including recombinanttechniques), cell biology, biochemistry, and immunology, which arewithin the skill of the art. Such conventional techniques includepolymer array synthesis, hybridization, ligation, and detection ofhybridization using a label. Specific illustrations of suitabletechniques can be had by reference to the example herein below. However,other equivalent conventional procedures can, of course, also be used.Such conventional techniques and descriptions can be found in standardlaboratory manuals such as Genome Analysis: A Laboratory Manual Series(Vols. I-IV), Using Antibodies: A Laboratory Manual, Cells: A LaboratoryManual, PCR Primer: A Laboratory Manual, and Molecular Cloning: ALaboratory Manual (all from Cold Spring Harbor Laboratory Press),Stryer, L. (1995) Biochemistry (4th Ed.) Freeman, New York, Gait,“Oligonucleotide Synthesis: A Practical Approach” 1984, IRL Press,London, Nelson and Cox (2000), Lehninger, Principles of Biochemistry3^(rd) Ed., W.H. Freeman Pub., New York, N.Y. and Berg et al. (2002)Biochemistry, 5^(th) Ed., W.H. Freeman Pub., New York, N.Y., all ofwhich are herein incorporated in their entirety by reference for allpurposes.

The present invention can employ solid substrates, including arrays insome preferred embodiments. Methods and techniques applicable to polymer(including protein) array synthesis have been described in U.S. Ser. No.09/536,841, WO 00/58516, U.S. Pat. Nos. 5,143,854, 5,242,974, 5,252,743,5,324,633, 5,384,261, 5,405,783, 5,424,186, 5,451,683, 5,482,867,5,491,074, 5,527,681, 5,550,215, 5,571,639, 5,578,832, 5,593,839,5,599,695, 5,624,711, 5,631,734, 5,795,716, 5,831,070, 5,837,832,5,856,101, 5,858,659, 5,936,324, 5,968,740, 5,974,164, 5,981,185,5,981,956, 6,025,601, 6,033,860, 6,040,193, 6,090,555, 6,136,269,6,269,846 and 6,428,752, in PCT Applications Nos. PCT/US99/00730(International Publication No. WO 99/36760) and PCT/US01/04285(International Publication No. WO 01/58593), which are all incorporatedherein by reference in their entirety for all purposes.

Patents that describe synthesis techniques in specific embodimentsinclude U.S. Pat. Nos. 5,412,087, 6,147,205, 6,262,216, 6,310,189,5,889,165, and 5,959,098. Nucleic acid arrays are described in many ofthe above patents, but the same techniques are applied to polypeptidearrays.

Nucleic acid arrays that are useful in the present invention includethose that are commercially available from Affymetrix (Santa Clara,Calif.) under the brand name GeneChip®. Example arrays are shown on thewebsite at affymetrix.com.

The present invention also contemplates many uses for polymers attachedto solid substrates. These uses include gene expression monitoring,profiling, library screening, genotyping and diagnostics. Geneexpression monitoring and profiling methods can be shown in U.S. Pat.Nos. 5,800,992, 6,013,449, 6,020,135, 6,033,860, 6,040,138, 6,177,248and 6,309,822. Genotyping and uses therefore are shown in U.S. Ser. Nos.10/442,021, 10/013,598 (U.S. Patent Application Publication20030036069), and U.S. Pat. Nos. 5,856,092, 6,300,063, 5,858,659,6,284,460, 6,361,947, 6,368,799 and 6,333,179. Other uses are embodiedin U.S. Pat. Nos. 5,871,928, 5,902,723, 6,045,996, 5,541,061, and6,197,506.

The present invention also contemplates sample preparation methods incertain preferred embodiments. Prior to or concurrent with genotyping,the genomic sample may be amplified by a variety of mechanisms, some ofwhich may employ PCR. See, for example, PCR Technology: Principles andApplications for DNA Amplification (Ed. H.A. Erlich, Freeman Press, NY,N.Y., 1992); PCR Protocols: A Guide to Methods and Applications (Eds.Innis, et al., Academic Press, San Diego, Calif., 1990); Mattila et al.,Nucleic Acids Res. 19, 4967 (1991); Eckert et al., PCR Methods andApplications 1, 17 (1991); PCR (Eds. McPherson et al., IRL Press,Oxford); and U.S. Pat. Nos. 4,683,202, 4,683,195, 4,800,159, 4,965,188,and 5,333,675, each of which is incorporated herein by reference intheir entireties for all purposes. The sample may be amplified on thearray. See, for example, U.S. Pat. No. 6,300,070 and U.S. Ser. No.09/513,300, which are incorporated herein by reference.

Other suitable amplification methods include the ligase chain reaction(LCR) (for example, Wu and Wallace, Genomics 4, 560 (1989), Landegren etal., Science 241, 1077 (1988) and Barringer et al. Gene 89:117 (1990)),transcription amplification (Kwoh et al., Proc. Natl. Acad. Sci. USA 86,1173 (1989) and WO88/10315), self-sustained sequence replication(Guatelli et al., Proc. Nat. Acad. Sci. USA, 87, 1874 (1990) andWO90/06995), selective amplification of target polynucleotide sequences(U.S. Pat. No. 6,410,276), consensus sequence primed polymerase chainreaction (CP-PCR) (U.S. Pat. No. 4,437,975), arbitrarily primedpolymerase chain reaction (AP-PCR) (U.S. Pat. Nos. 5,413,909, 5,861,245)and nucleic acid based sequence amplification (NABSA). (See, U.S. Pat.Nos. 5,409,818, 5,554,517, and 6,063,603, each of which is incorporatedherein by reference). Other amplification methods that may be used aredescribed in U.S. Pat. Nos. 5,242,794, 5,494,810, 4,988,617 and in U.S.Ser. No. 09/854,317, each of which is incorporated herein by reference.

Additional methods of sample preparation and techniques for reducing thecomplexity of a nucleic sample are described in Dong et al., GenomeResearch 11, 1418 (2001), in U.S. Pat. Nos. 6,361,947, 6,391,592 andU.S. Ser. Nos. 09/916,135, 09/920,491 (U.S. Patent ApplicationPublication 20030096235), Ser. No. 09/910,292 (U.S. Patent ApplicationPublication 20030082543), and Ser. No. 10/013,598.

Methods for conducting polynucleotide hybridization assays have beenwell developed in the art. Hybridization assay procedures and conditionswill vary depending on the application and are selected in accordancewith the general binding methods known including those referred to in:Maniatis et al. Molecular Cloning: A Laboratory Manual (2^(nd) Ed. ColdSpring Harbor, N.Y, 1989); Berger and Kimmel Methods in Enzymology, Vol.152, Guide to Molecular Cloning Techniques (Academic Press, Inc., SanDiego, Calif., 1987); Young and Davis, P.N.A.S, 80: 1194 (1983). Methodsand apparatus for carrying out repeated and controlled hybridizationreactions have been described in U.S. Pat. Nos. 5,871,928, 5,874,219,6,045,996 and 6,386,749, 6,391,623 each of which are incorporated hereinby reference

The present invention also contemplates signal detection ofhybridization between ligands in certain preferred embodiments. See U.S.Pat. Nos. 5,143,854, 5,578,832, 5,631,734, 5,834,758, 5,936,324,5,981,956, 6,025,601, 6,141,096, 6,185,030, 6,201,639, 6,218,803, and6,225,625 in U.S. Ser. No. 10/389,194 and in PCT ApplicationPCT/US99/06097 (published as WO99/47964), each of which also is herebyincorporated by reference in its entirety for all purposes.

Methods and apparatus for signal detection and processing of intensitydata are disclosed in, for example, U.S. Pat. Nos. 5,143,854, 5,547,839,5,578,832, 5,631,734, 5,800,992, 5,834,758, 5,856,092, 5,902,723,5,936,324, 5,981,956, 6,025,601, 6,090,555, 6,141,096, 6,185,030,6,201,639; 6,218,803; and 6,225,625, in U.S. Ser. Nos. 10/389,194,60/493,495 and in PCT Application PCT/US99/06097 (published asWO99/47964), each of which also is hereby incorporated by reference inits entirety for all purposes.

The practice of the present invention may also employ conventionalbiology methods, software and systems. Computer software products of theinvention typically include computer readable medium havingcomputer-executable instructions for performing the logic steps of themethod of the invention. Suitable computer readable medium includefloppy disk, CD-ROM/DVD/DVD-ROM, hard-disk drive, flash memory, ROM/RAM,magnetic tapes, etc. The computer-executable instructions may be writtenin a suitable computer language or combination of several languages.Basic computational biology methods are described in, for example,Setubal and Meidanis et al., Introduction to Computational BiologyMethods (PWS Publishing Company, Boston, 1997); Salzberg, Searles,Kasif, (Ed.), Computational Methods in Molecular Biology, (Elsevier,Amsterdam, 1998); Rashidi and Buehler, Bioinformatics Basics:Application in Biological Science and Medicine (CRC Press, London, 2000)and Ouelette and Bzevanis Bioinformatics: A Practical Guide for Analysisof Gene and Proteins (Wiley & Sons, Inc., 2^(nd) ed., 2001). See U.S.Pat. No. 6,420,108.

The present invention may also make use of various computer programproducts and software for a variety of purposes, such as probe design,management of data, analysis, and instrument operation. See, U.S. Pat.Nos. 5,593,839, 5,795,716, 5,733,729, 5,974,164, 6,066,454, 6,090,555,6,185,561, 6,188,783, 6,223,127, 6,229,911 and 6,308,170.

Additionally, the present invention may have preferred embodiments thatinclude methods for providing genetic information over networks such asthe Internet as shown in U.S. Ser. Nos. 10/197,621, 10/063,559 (UnitedStates Publication No. 20020183936), Ser. Nos. 10/065,856, 10/065,868,10/328,818, 10/328,872, 10/423,403, and 60/482,389 and U.S. PatentPublication No. 20040049354.

b) Definitions

The term “array” as used herein refers to an intentionally createdcollection of molecules which can be prepared either synthetically orbiosynthetically. The molecules in the array can be identical ordifferent from each other. The array can assume a variety of formats,for example, libraries of soluble molecules; libraries of compoundstethered to resin beads, silica chips, or other solid supports.

The term “array plate” as used herein refers to a body having aplurality of arrays in which each microarray is separated by a physicalbarrier resistant to the passage of liquids and forming an area orspace, referred to as a well, capable of containing liquids in contactwith the probe array.

The term “biomonomer” as used herein refers to a single unit ofbiopolymer, which can be linked with the same or other biomonomers toform a biopolymer (for example, a single amino acid or nucleotide withtwo linking groups one or both of which may have removable protectinggroups) or a single unit which is not part of a biopolymer. Thus, forexample, a nucleotide is a biomonomer within an oligonucleotidebiopolymer, and an amino acid is a biomonomer within a protein orpeptide biopolymer; avidin, biotin, antibodies, antibody fragments,etc., for example, are also biomonomers.

The term “biopolymer” or sometimes refer by “biological polymer” as usedherein is intended to mean repeating units of biological or chemicalmoieties. Representative biopolymers include, but are not limited to,nucleic acids, oligonucleotides, amino acids, proteins, peptides,hormones, oligosaccharides, lipids, glycolipids, lipopolysaccharides,phospholipids, synthetic analogues of the foregoing, including, but notlimited to, inverted nucleotides, peptide nucleic acids, Meta-DNA, andcombinations of the above.

The term “biopolymer synthesis” as used herein is intended to encompassthe synthetic production, both organic and inorganic, of a biopolymer.Related to a bioploymer is a “biomonomer”.

The term “cartridge” as used herein refers to a body forming an area orspace referred to as a well wherein a microarray is contained andseparated from the passage of liquids.

The term “clamping plate” as used herein refers to a device used forfastening two or more parts.

The term “combinatorial synthesis strategy” as used herein refers to anordered strategy for parallel synthesis of diverse polymer sequences bysequential addition of reagents which may be represented by a reactantmatrix and a switch matrix, the product of which is a product matrix. Areactant matrix is al column by m row matrix of the building blocks tobe added. The switch matrix is all or a subset of the binary numbers,preferably ordered, between 1 and m arranged in columns. A “binarystrategy” is one in which at least two successive steps illuminate aportion, often half, of a region of interest on the substrate. In abinary synthesis strategy, all possible compounds which can be formedfrom an ordered set of reactants are formed. In most preferredembodiments, binary synthesis refers to a synthesis strategy which alsofactors a previous addition step. For example, a strategy in which aswitch matrix for a masking strategy halves regions that were previouslyilluminated, and then illuminating about half of the previouslyilluminated region and protecting the remaining half (while alsoprotecting about half of previously protected regions and illuminatingabout half of previously protected regions). It will be recognized thatbinary rounds may be interspersed with non-binary rounds and that only aportion of a substrate may be subjected to a binary scheme. Acombinatorial “masking” strategy is a synthesis which uses light orother spatially selective deprotecting or activating agents to removeprotecting groups from materials for addition of other materials such asamino acids.

The term “complementary” as used herein refers to the hybridization orbase pairing between nucleotides or nucleic acids, such as, forinstance, between the two strands of a double stranded DNA molecule orbetween an oligonucleotide primer and a primer binding site on a singlestranded nucleic acid to be sequenced or amplified. Complementarynucleotides are, generally, A and T (or A and U), or C and G. Two singlestranded RNA or DNA molecules are said to be complementary when thenucleotides of one strand, optimally aligned and compared and withappropriate nucleotide insertions or deletions, pair with at least about80% of the nucleotides of the other strand, usually at least about 90%to 95%, and more preferably from about 98 to 100%. Alternatively,complementarity exists when an RNA or DNA strand will hybridize underselective hybridization conditions to its complement. Typically,selective hybridization will occur when there is at least about 65%complementary over a stretch of at least 14 to 25 nucleotides,preferably at least about 75%, more preferably at least about 90%complementary. See, M. Kanehisa Nucleic Acids Res. 12:203 (1984),incorporated herein by reference.

The term “effective amount” as used herein refers to an amountsufficient to induce a desired result.

The term “excitation energy” as used herein refers to energy used toenergize a detectable label for detection, for example illuminating afluorescent label. Devices for this use include coherent light or noncoherent light, such as lasers, UV light, light emitting diodes, anincandescent light source, or any other light or other electromagneticsource of energy having a wavelength in the excitation band of anexcitable label, or capable of providing detectable transmitted,reflective, or diffused radiation.

The term “gaskets or o-ring” as used herein refers to any of a widevariety of seals or packings used between joined parts to prevent theescape of a gas or fluid. Gaskets or o-rings can be made of materialssuch as elastomer.

The term “genome” as used herein is all the genetic material in thechromosomes of an organism. DNA derived from the genetic material in thechromosomes of a particular organism is genomic DNA. A genomic libraryis a collection of clones made from a set of randomly generatedoverlapping DNA fragments representing the entire genome of an organism.

The term “hybridization” as used herein refers to the process in whichtwo single- stranded polynucleotides bind noncovalently to form a stabledouble-stranded polynucleotide; triple-stranded hybridization is alsotheoretically possible. The resulting (usually) double-strandedpolynucleotide is a “hybrid.” The proportion of the population ofpolynucleotides that forms stable hybrids is referred to herein as the“degree of hybridization.” Hybridizations are usually performed understringent conditions, for example, at a salt concentration of no morethan 1 M and a temperature of at least 25° C. For example, conditions of5× SSPE (750 mM NaCl, 50 mM NaPhosphate, 5 mM EDTA, pH 7.4) and atemperature of 25-30° C. are suitable for allele-specific probehybridizations. For stringent conditions, see, for example, Sambrook,Fritsche and Maniatis. “Molecular Cloning A laboratory Manual” 2^(nd)Ed. Cold Spring Harbor Press (1989) which is hereby incorporated byreference in its entirety for all purposes above.

The term “hybridization conditions” as used herein will typicallyinclude salt concentrations of less than about 1M, more usually lessthan about 500 mM and preferably less than about 200 mM. Hybridizationtemperatures can be as low as 5° C., but are typically greater than 22°C., more typically greater than about 30° C., and preferably in excessof about 37° C. Longer fragments may require higher hybridizationtemperatures for specific hybridization. As other factors may affect thestringency of hybridization, including base composition and length ofthe complementary strands, presence of organic solvents and extent ofbase mismatching, the combination of parameters is more important thanthe absolute measure of any one alone.

The term “hybridization probes” as used herein are oligonucleotidescapable of binding in a base-specific manner to a complementary strandof nucleic acid. Such probes include peptide nucleic acids, as describedin Nielsen et al., Science 254, 1497-1500 (1991), and other nucleic acidanalogs and nucleic acid mimetics.

The term “hybridizing specifically to” as used herein refers to thebinding, duplexing, or hybridizing of a molecule only to a particularnucleotide sequence or sequences under stringent conditions when thatsequence is present in a complex mixture (for example, total cellular)DNA or RNA.

The term “isolated nucleic acid” as used herein means an object speciesinvention that is the predominant species present (i.e., on a molarbasis it is more abundant than any other individual species in thecomposition). Preferably, an isolated nucleic acid comprises at leastabout 50, 80 or 90% (on a molar basis) of all macromolecular speciespresent. Most preferably, the object species is purified to essentialhomogeneity (contaminant species cannot be detected in the compositionby conventional detection methods).

The term “label” as used herein refers to a luminescent label, a lightscattering label or a radioactive label. Fluorescent labels include,inter alia, the commercially available fluorescein phosphoramidites suchas Fluoreprime (Pharmacia), Fluoredite (Millipore) and FAM (ABI). SeeU.S. Pat. No. 6,287,778.

“Labeled moiety” refers to a moiety capable of being detected by thevarious methods discussed herein or known in the art.

The term “ligand” as used herein refers to a molecule that is recognizedby a particular receptor. The agent bound by or reacting with a receptoris called a “ligand,” a term which is definitionally meaningful only interms of its counterpart receptor. The term “ligand” does not imply anyparticular molecular size or other structural or compositional featureother than that the substance in question is capable of binding orotherwise interacting with the receptor. Also, a ligand may serve eitheras the natural ligand to which the receptor binds, or as a functionalanalogue that may act as an agonist or antagonist. Examples of ligandsthat can be investigated by this invention include, but are notrestricted to, agonists and antagonists for cell membrane receptors,toxins and venoms, viral epitopes, hormones (for example, opiates,steroids, etc.), hormone receptors, peptides, enzymes, enzymesubstrates, substrate analogs, transition state analogs, cofactors,drugs, proteins, and antibodies.

The term “microtiter plates” as used herein refers to arrays of discretewells that come in standard formats (96, 384 and 1536 wells) which areused for examination of the physical, chemical or biologicalcharacteristics of a quantity of samples in parallel.

The term “mixed population” or sometimes refer by “complex population”as used herein refers to any sample containing both desired andundesired nucleic acids. As a non-limiting example, a complex populationof nucleic acids may be total genomic DNA, total genomic RNA or acombination thereof. Moreover, a complex population of nucleic acids mayhave been enriched for a given population but include other undesirablepopulations. For example, a complex population of nucleic acids may be asample which has been enriched for desired messenger RNA (mRNA)sequences but still includes some undesired ribosomal RNA sequences(rRNA).

The term “mRNA” or sometimes refer by “mRNA transcripts” as used herein,include, but not limited to pre-mRNA transcript(s), transcriptprocessing intermediates, mature mRNA(s) ready for translation andtranscripts of the gene or genes, or nucleic acids derived from the mRNAtranscript(s). Transcript processing may include splicing, editing anddegradation. As used herein, a nucleic acid derived from an mRNAtranscript refers to a nucleic acid for whose synthesis the mRNAtranscript or a subsequence thereof has ultimately served as a template.Thus, a cDNA reverse transcribed from an mRNA, an RNA transcribed fromthat cDNA, a DNA amplified from the cDNA, an RNA transcribed from theamplified DNA, etc., are all derived from the mRNA transcript anddetection of such derived products is indicative of the presence and/orabundance of the original transcript in a sample. Thus, mRNA derivedsamples include, but are not limited to, mRNA transcripts of the gene orgenes, cDNA reverse transcribed from the mRNA, cRNA transcribed from thecDNA, DNA amplified from the genes, RNA transcribed from amplified DNA,and the like.

The term “nucleic acid library” or sometimes refer by “array” as usedherein refers to an intentionally created collection of nucleic acidswhich can be prepared either synthetically or biosynthetically andscreened for biological activity in a variety of different formats (forexample, libraries of soluble molecules; and libraries of oligostethered to resin beads, silica chips, or other solid supports).Additionally, the term “array” is meant to include those libraries ofnucleic acids which can be prepared by spotting nucleic acids ofessentially any length (for example, from 1 to about 1000 nucleotidemonomers in length) onto a substrate. The term “nucleic acid” as usedherein refers to a polymeric form of nucleotides of any length, eitherribonucleotides, deoxyribonucleotides or peptide nucleic acids (PNAs),that comprise purine and pyrimidine bases, or other natural, chemicallyor biochemically modified, non-natural, or derivatized nucleotide bases.The backbone of the polynucleotide can comprise sugars and phosphategroups, as may typically be found in RNA or DNA, or modified orsubstituted sugar or phosphate groups. A polynucleotide may comprisemodified nucleotides, such as methylated nucleotides and nucleotideanalogs. The sequence of nucleotides may be interrupted bynon-nucleotide components. Thus the terms nucleoside, nucleotide,deoxynucleoside and deoxynucleotide generally include analogs such asthose described herein. These analogs are those molecules having somestructural features in common with a naturally occurring nucleoside ornucleotide such that when incorporated into a nucleic acid oroligonucleoside sequence, they allow hybridization with a naturallyoccurring nucleic acid sequence in solution. Typically, these analogsare derived from naturally occurring nucleosides and nucleotides byreplacing and/or modifying the base, the ribose or the phosphodiestermoiety. The changes can be tailor made to stabilize or destabilizehybrid formation or enhance the specificity of hybridization with acomplementary nucleic acid sequence as desired.

The term “nucleic acids” as used herein may include any polymer oroligomer of pyrimidine and purine bases, preferably cytosine, thymine,and uracil, and adenine and guanine, respectively. See Albert L.Lehninger, PRINCIPLES OF BIOCHEMISTRY, at 793-800 (Worth Pub. 1982).Indeed, the present invention contemplates any deoxyribonucleotide,ribonucleotide or peptide nucleic acid component, and any chemicalvariants thereof, such as methylated, hydroxymethylated or glucosylatedforms of these bases, and the like. The polymers or oligomers may beheterogeneous or homogeneous in composition, and may be isolated fromnaturally-occurring sources or may be artificially or syntheticallyproduced. In addition, the nucleic acids may be DNA or RNA, or a mixturethereof, and may exist permanently or transitionally in single-strandedor double-stranded form, including homoduplex, heteroduplex, and hybridstates.

The term “oligonucleotide” or sometimes refer by “polynucleotide” asused herein refers to a nucleic acid ranging from at least 2, preferableat least 8, and more preferably at least 20 nucleotides in length or acompound that specifically hybridizes to a polynucleotide.Polynucleotides of the present invention include sequences ofdeoxyribonucleic acid (DNA) or ribonucleic acid (RNA) which may beisolated from natural sources, recombinantly produced or artificiallysynthesized and mimetics thereof. A further example of a polynucleotideof the present invention may be peptide nucleic acid (PNA). Theinvention also encompasses situations in which there is a nontraditionalbase pairing such as Hoogsteen base pairing which has been identified incertain tRNA molecules and postulated to exist in a triple helix.“Polynucleotide” and “oligonucleotide” are used interchangeably in thisapplication.

The term “polymorphism” as used herein refers to the occurrence of twoor more genetically determined alternative sequences or alleles in apopulation. A polymorphic marker or site is the locus at whichdivergence occurs. Preferred markers have at least two alleles, eachoccurring at frequency of greater than 1%, and more preferably greaterthan 10% or 20% of a selected population. A polymorphism may compriseone or more base changes, an insertion, a repeat, or a deletion. Apolymorphic locus may be as small as one base pair. Polymorphic markersinclude restriction fragment length polymorphisms, variable number oftandem repeats (VNTR's), hypervariable regions, mini satellites,dinucleotide repeats, trinucleotide repeats, tetranucleotide repeats,simple sequence repeats, and insertion elements such as Alu. The firstidentified allelic form is arbitrarily designated as the reference formand other allelic forms are designated as alternative or variantalleles. The allelic form occurring most frequently in a selectedpopulation is sometimes referred to as the wildtype form. Diploidorganisms may be homozygous or heterozygous for allelic forms. Adiallelic polymorphism has two forms. A triallelic polymorphism hasthree forms. Single nucleotide polymorphisms (SNPs) are included inpolymorphisms.

The term “primer” as used herein refers to a single-strandedoligonucleotide capable of acting as a point of initiation fortemplate-directed DNA synthesis under suitable conditions for example,buffer and temperature, in the presence of four different nucleosidetriphosphates and an agent for polymerization, such as, for example, DNAor RNA polymerase or reverse transcriptase. The length of the primer, inany given case, depends on, for example, the intended use of the primer,and generally ranges from 15 to 30 nucleotides. Short primer moleculesgenerally require cooler temperatures to form sufficiently stable hybridcomplexes with the template. A primer need not reflect the exactsequence of the template but must be sufficiently complementary tohybridize with such template. The primer site is the area of thetemplate to which a primer hybridizes. The primer pair is a set ofprimers including a 5′ upstream primer that hybridizes with the 5′ endof the sequence to be amplified and a 3′ downstream primer thathybridizes with the complement of the 3′ end of the sequence to beamplified.

The term “probe” as used herein refers to a surface-immobilized moleculethat can be recognized by a particular target. See U.S. Pat. No.6,582,908 for an example of arrays having all possible combinations ofprobes with 10, 12, and more bases. Examples of probes that can beinvestigated by this invention include, but are not restricted to,agonists and antagonists for cell membrane receptors, toxins and venoms,viral epitopes, hormones (for example, opioid peptides, steroids, etc.),hormone receptors, peptides, enzymes, enzyme substrates, cofactors,drugs, lectins, sugars, oligonucleotides, nucleic acids,oligosaccharides, proteins, and monoclonal antibodies.

The term “reader” or “plate reader” as used herein refers to a devicewhich is used to identify hybridization events on an array, such as thehybridization between a nucleic acid probe on the array and afluorescently labeled target. Readers are known in the art and arecommercially available through Affymetrix, Santa Clara Calif. and othercompanies. Generally, they involve the use of an excitation energy (suchas a laser) to illuminate a fluorescently labeled target nucleic acidthat has hybridized to the probe. Then, the reemitted radiation (at adifferent wavelength than the excitation energy) is detected usingdevices such as a CCD, PMT, photodiode, or similar devices to registerthe collected emissions. See U.S. Pat. No. 6,225,625.

The term “receptor” as used herein refers to a molecule that has anaffinity for a given ligand. Receptors may be naturally-occurring ormanmade molecules. Also, they can be employed in their unaltered stateor as aggregates with other species. Receptors may be attached,covalently or noncovalently, to a binding member, either directly or viaa specific binding substance. Examples of receptors which can beemployed by this invention include, but are not restricted to,antibodies, cell membrane receptors, monoclonal antibodies and antiserareactive with specific antigenic determinants (such as on viruses, cellsor other materials), drugs, polynucleotides, nucleic acids, peptides,cofactors, lectins, sugars, polysaccharides, cells, cellular membranes,and organelles. Receptors are sometimes referred to in the art asanti-ligands. As the term receptors is used herein, no difference inmeaning is intended. A “Ligand Receptor Pair” is formed when twomacromolecules have combined through molecular recognition to form acomplex. Other examples of receptors which can be investigated by thisinvention include but are not restricted to those molecules shown inU.S. Pat. No. 5,143,854, which is hereby incorporated by reference inits entirety.

The term “solid support”, “support”, and “substrate” as used herein areused interchangeably and refer to a material or group of materialshaving a rigid or semi-rigid surface or surfaces. In many embodiments,at least one surface of the solid support will be substantially flat,although in some embodiments it may be desirable to physically separatesynthesis regions for different compounds with, for example, wells,raised regions, pins, etched trenches, or the like. According to otherembodiments, the solid support(s) will take the form of beads, resins,gels, microspheres, or other geometric configurations. See U.S. Pat. No.5,744,305 for exemplary substrates.

The term “surface” or “active probe surface” or “target surface” as usedherein refers to the area of the microarray to be analyzed withreagents.

The term “target” as used herein refers to a molecule that has anaffinity for a given probe. Targets may be naturally-occurring orman-made molecules. Also, they can be employed in their unaltered stateor as aggregates with other species. Targets may be attached, covalentlyor noncovalently, to a binding member, either directly or via a specificbinding substance. Examples of targets which can be employed by thisinvention include, but are not restricted to, antibodies, cell membranereceptors, monoclonal antibodies and antisera reactive with specificantigenic determinants (such as on viruses, cells or other materials),drugs, oligonucleotides, nucleic acids, peptides, cofactors, lectins,sugars, polysaccharides, cells, cellular membranes, and organelles.Targets are sometimes referred to in the art as anti-probes. As the termtargets is used herein, no difference in meaning is intended. A “ProbeTarget Pair” is formed when two macromolecules have combined throughmolecular recognition to form a complex.

The term “variance” as used herein generally refers to a value that is ameasure of the dispersion of data. For example, it will be appreciatedby those skilled in the relevant art that, variance may be defined asthe mean of the square of the differences between the samples and can bemathematically represented as:

$\sigma^{2} = \frac{{\Sigma ( {X - \overset{\_}{X}} )}^{2}}{n - 1}$

where, X is equal to a particular value that could for instance be anemission intensity value for a probe feature, X is equal to the mean ofall the values and n is equal to the total number of values. The mean isequal to the sum of all the observation values divided by the number ofobservations.

The term “wafer” as used herein refers to a substrate having surface towhich a plurality of arrays are bound. In a preferred embodiment, thearrays are synthesized on the surface of the substrate to createmultiple arrays that are physically separate. In one preferredembodiment of a wafer, the arrays are physically separated by a distanceof at least about 0.1, 0.25, 0.5, 1 or 1.5 millimeters. The arrays thatare on the wafer may be identical, each one may be different, or theremay be some combination thereof. Particularly preferred wafers are about8″×8″ and are made using the photolithographic process.

C. Analysis of Gene Expression on Exon Arrays I. Alternative TranscriptAnalysis Methods for Exon Arrays

With exon arrays like the GeneChip® Human Exon 1.0 ST Array, researcherscan examine the transcriptional profile of an entire gene. Being able togather data for each individual exon enables the investigation ofphenomena such as alternative splicing, alternative promoter usage, andalternative termination while also providing more probe level data todetermine the overall rate of expression a particular locus. Methods ofdetecting relative changes in alternative transcript forms are disclosedherein.

An ANOVA based method based on a Splicing Index equal to the ratio ofexon signal to gene signal works well in two test data sets: one atissue panel and the other a set of cancer and normal from the sametissue. For the tissue panel data set, a correlation based method,Robust PAC, performs well. Other methods for alternative splicinganalysis bases on ANOVA are disclosed in U.S. patent application Ser.No. 11/069,720.

One of the most well studied alternative splicing events is thealternative utilization of “cassette exons”. The performance ofdifferent mathematical methods designed to detect alternatively splicedcassette exons from the exon array data was evaluated. First, themethod's robustness was tested in an artificial data set by constructing“genes” comprised of a group of probe selection regions that are bestsuggested by annotation evidence as biologically real and are alwaysjointly expressed in a sample data set. Each such group of probeselection regions is referred to as “core constitutive exons” of that“gene.” The alternative splice events were simulated in each such “gene”by substituting exons from alternative regions of the genome andconstructed ROC curves as a way to measure the performance of thedifferent methods.

The data set included a set of core constitutive exons, a tissue paneland a colon cancer panel. To generate a set of exons that appear to beconstitutively included in all transcripts (core constitutive exons)splicing graphs were generated as described (Sugnet, et al, Pac SympBiocomput., p. 66-77, 2004) using all human RefSeq, mRNA and ESTsequences. The graphs were pruned by requiring that each exon be eitheralso present in mouse cDNAs or present multiple times in human cDNAs.The graphs were then searched to identify exons that were constitutivelyincluded in the pruned graph. To ensure that an exon is likely to bepresent in all transcripts a minimum number of 5 mRNAs or 10 ESTs (or acombination of the two) had to contain the exon. Exons making up thecore constitutive genes tend to be well-expressed; their median signalis greater than 10 times the median signal of background probesets inboth data sets.

The tissue panel data set consists of 11 human tissue samples. Thehybridization cocktail for each sample was scanned in triplicate, usingthree arrays per sample for a total of 33 scans. In this study wequantile normalized the triplicates together, then multiplicativelyscaled the resulting scans so that they all had the same median (“medianscaling”). All features, whether genomic or non-genomic are included inboth normalizations. The alternative splice set was generated by takingthe probes of one exon in each gene (the first in each set) andsubstituting them into another gene. The probes of the real exon in thesubstituted gene were moved to yet another gene, and after repeatingover all 5,800 genes, each gene had an “exon” that did not belong to it.

This approach works as variation within replicate groups is not due tobiological variation and hence, if the core constitutive genes havedetectable expression in one or more of the tissues, they will likelyhave different expression in the different tissues. Hence thealternative splice set will each contain an exon that will behavedifferently across tissues. The major caveat in creating thisalternative splice set is that it is not clear how well this proceduremimics true biological alternative splicing behavior.

The Colon Cancer data set consists of 18 biologically distinct samplesarranged in 9 pairs. Each pair is a normal colon tissue sample and acolon cancer sample from each of 9 different individuals. There are noreplicate scans. All scans were normalized using median scaling. For thecolon cancer set, biological variation in signal is large enough withinnormal and tumor groups so that around 98% of the genes did not showsignificant differences in means of signal between the group of cancerand normal samples (p-value>0.05). Accordingly, the approach in creatingan Alternative splice set as in the tissue panel data will not workwell; replacing an exon by an exon from another gene will not inducedetectable signal changes from one group to the other and hence methodsof alternative splice detection based on changes in signal of exonnormalized by gene will not detect anything. Instead, we took adifferent approach: a “background” set of probesets was generated byrandomly selecting 5,800 probesets from probesets with mean signalacross all 18 samples less than the 45^(th) percentile of mean signal ofeach probeset across all 18 samples.

An alternative splice data set was simulated by replacing probeintensities in the first exon in each gene with probe intensities from a“background” probeset for the normal scans only. Tumor scans were leftuntouched. We matched the number of probes in the replacement exon, asPLIER (PLIER: An M-Estimator for Expression Array, Hubbell, E.Affymetrix White Paper, 2005) signal is calculated across all samples byusing probe intensities. Replacement probes need not have the same GCcontent, and replacement probe intensities in each sample were adjustedby the ratio of surrogate intensity mismatch of original to replacementbased on GC content.

Exons making up the core constitutive genes have much higher overallsignal than the background probesets on the chip: the median signal(2.5) of background probesets is around the 20^(th) percentile of allprobeset signals and the median signal (35) of the first exon is aroundthe 70^(th) percentile of all probeset signals.

Several different methods for splice detection are disclosed:Pattern-Based Correlation (PAC) and two ANOVA methods: MicroarrayDetection of Alternative Splicing (MIDAS) presented here, and ANOSVA(Cline et al, Bioinformatics 21(suppl 1):i107-i115, 2005). Theirsuitability to any particular data set will depend on the structure ofthe data set and the goals of the experiment. Data set explorationpreferably uses at least two methods, including variations on thedisclosed methods.

In general, preferred splice detection methods include the followingfeatures. Under the null hypothesis, exons or probes comprising a geneare assumed to be proportional to each other across different samples. Amodel is fit that predicts probe or exon response under the nullhypothesis. A statistic is constructed that measures how deviant thedata is from the model. This statistic is used to construct a p-value.

Splicing Index

A Splicing Index captures the basic metric for the analysis ofalternative splicing. Specifically it is a measure of how much exonspecific expression (with gene induction factored out) differs betweentwo samples. The first step is to normalize the exon level signals tothe gene level signals.

$\begin{matrix}{n_{i,j,k} = \frac{e_{i,j,k}}{g_{j,k}}} & {{Equation}\mspace{14mu} 1}\end{matrix}$

where e_(i,j,k) is the exon signal estimate of the i exon, j experiment,and k gene. g_(i,k) is the gene level signal estimate of the jexperiment and k gene. The splicing index is then the ratio ofnormalized exon signal estimates from one sample or set of samplesrelative to another. For example, in the colon cancer data set asplicing index could be established for each gene by taking the mediann_(i,j,k) for the cancer samples and dividing by the median n_(i,j,k)for the normal samples. Alternatively one might want to calculate asplicing ratio for each paired normal/cancer sample and then report themedian ratio. A similar approach is described in Clark et al. Science296:907-10 (2002).

The PAC approach assumes that in the absence of splicing, exonexpression follows gene expressions across samples using the followingmodel:

e_(i,j,k)=α_(i,k)g_(j,k)   Equation 2

where e_(i,j,k) is the signal of the i-th exon of the j-th sample of thek-th gene, g_(j,k) is the signal of k-th gene in the j-th sample, andα_(i,k) is the ratio of exon i signal to its gene signal.

Preferably use a robust measure of gene signal and correlate signal ofeach exon with this signal. Low correlations are indications ofalternative splicing. The robust measure of gene signal allows multipleexons to not track the overall gene, so long as they remain in a “smallenough” minority. In preferred aspects PLIER is used.

An important class of experiments asks the following question: Is therean alternative splice variant present in one group out of two groups ofsamples. PAC has the problem that it will fail in two-sample cases, ascorrelation will always be +1 (exon response agrees in direction fromone sample type to the other) or −1 (exon response disagrees indirection); this is expected to happen no matter how small the changeactually is; in the rare event that there is no numerical changecorrelation would be zero. Hence PAC is better suited to experimentswith a relatively large number of sample types.

MIDAS

In another aspect an alternative ANOVA based method based on measuringdifferences between exon level signal and aggregate gene level signalmay be used. This method is referred to herein as Microarray Detectionof Alternative Splicing or MIDAS. In general the underlying concept ofMIDAS is as follows: (1) PLIER is used to generate a robust estimate ofgene-level signal by using data from all features in all exons in thegene. This signal has the virtue that it will be robust against exonsthat exhibit anomalous signal across samples, whether they benon-expressing probe sets, probe sets that are incorrectly assigned to agene, or are true alternative splice exons. Since this array has nomismatch probes a surrogate estimate of mismatch is generated and usedby using the median intensity of all antigenomic probes on the samearray that have the same GC content as the probe. (2) PLIER is used tosimilarly generate an estimate of signal for each exon. (3) Under thenull hypothesis of no alternative splicing for an exon, the differencebetween the logged signal for the exon and its gene is expected to be aconstant across all samples. In other words, it is expected thatobserved signal from each exon will have a constant ratio with observedsignal from its gene.

A detection metric or statistic will preferably be based on log ratio ofexon signal to gene signal, i.e, the difference between logged signal ofeach exon and its gene. In preferred aspects a small constant is addedbefore logging to stabilize variance.

MIDAS Relation to ANOVA

Statistics based on the Splicing Index can be calculated either for eachgene or for each exon. A model for possible splicing is:

e_(i,j,k)=α_(i,k)p_(i,j,k)g_(j,k)   Equation 3

where e_(i,j,k) is the signal of the i-th exon of the j-th sample of thek-th gene, g_(j,k) is the signal of k-th gene in the j-th sample,α_(i,k) is the ratio of exon i signal to its gene signal in the samplewhere it is maximally expressed, and 0≤p_(i,j,k)≤1 is a Splicing Indexthat estimates the proportionate expression of this exon of this gene intissue j. Dividing both sides by g_(j,k) to factor out effects due togene induction and taking logs reduces this to an additive model(ignoring the possibility of zero signal):

log(e _(i,j,k) /g _(j,k))=log(e _(i,j,k))−log(g_(j,k))=log(α_(i,k))+log(p _(i,j,k))   Equation 4

Exon-Level Detection vs. Gene Level DetectionGene-level MIDAS is a 2-way ANOVA that includes an error term andpossible interactions comparing:

log(e _(i,j,k))−log(g _(j,k))=log(α_(i,k))+log(p_(i,j,k))+γ_(i,j,k)+ε_(i,j,k)   Equation 5

to the reduced model:

log(e _(i,j,k))−log(g _(j,k))=log(α_(i,k))+ε_(i,j,k)   Equation 6

So we ask the question: are effects other than exon effects present;i.e., test log(p_(i,j,k))=γ_(i,j,k)=0 across samples and exons.Exon-level MIDAS considers the situation an exon at a time. In thisnarrow view log(α_(i,k)) is constant and hence it is appropriate to useclassical 1-way ANOVA that compares the full model:

log(e _(i,j,k))−log(g _(j,k))=log(α_(i,k))+log(p _(i,j,k))+ε_(i,j,k)  Equation 7

to:

log(e _(i,j,k))−log(g _(j,k))=log(α_(i,k))+ε_(i,j,k)   Equation 8

to test the hypothesis of no alternative splicing by testing for theconstant effects model log(p_(i,j,k))=0 for all k samples.

As discussed above, the log model is inappropriate when exons or geneare not expressed. In practice, stabilizing variance by adding aconstant (as we do here and further discussed below) helps to reduce thefalse positive rate. Both PAC and MIDAS show considerable improvement inthe ROC curves when using exon-level detection over gene-leveldetection.

ANOSVA

A detailed description of the ANOSVA can be found in (Cline et al,2005). In short, ANOSVA assumes all probe responses can be modeled asproportional to each other across different samples in the absence ofalternative splicing; this reduces to an additive ANOVA model aftertaking logs. The ANOSVA method then tests for an alternative hypothesisof non-zero interactions between samples and exons using an F-statistic.Preliminary evaluation of ANOSVA on exon array data did not yield goodperformance for exon array data.

DECONV

Another algorithm for estimating relative concentrations of differentsplice variants is described in (Wang, et al, Bioinformatics 200319(suppl 1):i315-i322), where it is applied to data on an experimentalmicroarray. If we define a set of splice variants each consisting of allexons but one, this algorithm can be recast into a model for probe-levelintensities x_(h,i,j,k) similar in form to Equation 2 above, but with anextra multiplicative term denoting the probe affinity a_(h,i,k) of probeh of exon i of gene k in tissue j:

x _(h,i,j,k) =a _(h,i,k)α_(i,k) p _(i,j,k) g _(j,k)+ε_(h,i,j,k)  Equation 9

Parameters are estimated using an iterative maximum likelihoodestimation method, including p_(i,j,k). For any one gene k, the relativeratios of the p_(i,j,k) across tissues j represent relativeconcentration estimates of exons. Alternative splicing would be presentif these relative ratios deviate from each other in a statisticallysignificant way. This model differs from MIDAS in that error is treatedadditively rather than multiplicatively, and gene signal is estimateddirectly from the MLE rather than from a separate PLIER fit.

Non-Expressing Probesets and Genes

Any Splicing Index based on the ration of exon signal to gene signalfinds probesets whose response data disagrees with a model wheresplicing is assumed absent. However, this can occur in cases with noalternative splicing: In an exploratory array like the Human Exon Array,many probesets are based on gene prediction algorithms and may beinterrogating regions that are not actually transcribed in a particularsample, or may be interrogating mistaken genomic assemblies, and hencedo not correspond to biological reality. Such probesets will typicallyhave low probe intensities corresponding to noise and the exon signalill not generally track the gene signal. When the gene and probeset docorrespond to biological reality, the gene may be poorly expressed orunexpressed in the biological samples. Both the gene and the probe setwill have low probe intensities corresponding to noise and the ratio ofexon signal to gene signal will have no meaning. Hence it is importantthat any statistical method based on the ratio of exon signal to genesignal should be able to handle low intensity probe sets.

Handling of probesets or genes not expressed in sample is an active areaof research and while there is no specific assessment of the methods inthe context of the model breakdown above, we took the following steps toguard against model breakdown:

Stabilize Variance

Most of the models discussed use the log ratio of exon signal to genesignal and for these we stabilized variance by adding a constant to thePLIER signal estimate before taking the ratio. This trades bias forvariance. We chose a constant approximately equal to the 20%-ile ofprobeset signals. This appears to be well within the background level ofprobeset signals and since the core constitutive genes tend to bewell-expressed, the bias will be relatively small.

Use Exon-Level Models

Probe sets that show only noise against a backdrop of true variation ingene signal will tend to generate higher values using MIDASF-statistics. If the detection statistic is calculated for the gene as awhole, then such exons must be removed before this calculation otherwisethe gene will be flagged as a candidate for alternative splice events.The robust nature of PLIER may be exploited to obtain estimates of genesignal to automatically filter out such probesets from the gene signal;this favors an exon-level approach.

Performance Evaluation

A ROC curve measures how well a statistic differentiates truealternatives from false positives. To do this exactly requires a knownset that do not exhibit alternative splicing (the null set) to becompared with a known set that do exhibit alternative splicing (thealternative set). However, in our situation, we do not have completeknowledge, and hence the null set is likely contaminated by some realalternative splicing.

The effect of mixing some true alternative splice data into the null setwill tend to cause the ROC curves to be constrained to a diamond(Bourgon, in preparation) around the diagonal of ROC plot and hencereduce the differences between statistics for large p-values. Mixingnull data into the alternative set also affects the shape of thediamond; in this case the differences between statistics for smallp-values will be affected.

It is expected that small p-values are more important than largep-values in assessing performance of different statistics that detectalternative splicing, and hence the error in construction of the nullset is unimportant in our conclusion. In any case, qualitatively theempirical ROC curves are well-behaved; errors in construction of thenull set probably do not affect p-values less than 0.1

Tissue Data Set Performance

Results vary considerably depending on the data set. With no filteringfor the tissue panel both PAC and MIDAS perform equivalently.

Colon Cancer/Normal Data Set Performance

On the colon cancer data set, PAC is not applicable, and MIDAS reducesto a 2-sample t-test. (See FIG. 2. ROC curve for Colon samples (Normalvs. Tumor) using MIDAS).

Methods Stratified by Number of Exons per Gene

A histogram of the distribution of exons is shown below in FIG. 3:

We split exons into approximately 5 equal groups binned by number ofexons per gene

(Error! Reference Source Not Found.):

Min exons/gene  3+  5+  10+ Max exons/gene 5 10 15 Number genes 1762  2252  1786 We see fairly similar ROC curves, with discrimination improving withlarger number of exons per gene, possibly because the PLIER estimate ofgene signal is more stable with a larger number of exons.

In practice, using exon-level detection instead of gene-level detectionforces the issue of multiple comparisons, where the larger the number ofexons in the gene, the greater the chance of high statisticalsignificance by chance alone.

In practice it is also desirable to use p-values. Using MIDAS on thetissue panel data, the ANOVA F-statistic gives p-values that are about 3times too low (FIG. 4) Comparison of Empirical p-Values with Theoreticalp-Values), so in practice one would not want to take p-values tooseriously. Different structure in replicates and sampling might not givethe same relation.

Using the method on RefSeq genes, highly significant p-values can beassociated with exons and genes that consistently exhibit low expressionacross sample types (model failures as discussed below). P-values canalso be high when an exon from one gene is incorrectly assigned toanother gene. While the ratio of exon signal to gene is variable,observing a very high ratio in one or more of the sample types is anindication that an error of this type may have occurred.

Exon Array Background Correction: Similar performance is expected indetection of differential changes for most applications on exon arraysusing the background probe collection (BGP) as opposed to specificmismatch probes (MM). (For example, the GeneChip® Human Exon 1.0 STArray has a “genomic” and an “antigenomic” background probe collection.See below for more information on the background probe collections.)

The use of BGP probes to perform a PM-GCBG correction is describedherein. PM-GCBG refers to the use of the median BGP intensity for BGPprobes with the same GC content as the perfect match (PM) probe.

MM probes provide controls similar in sequence and location which areuseful in removing bias; this takes one probe per informative probe.Because the MM probes hybridize to the perfect match (PM) target, whenMM intensities are subtracted from PM intensities the PM response totrue target slightly reduced. BGP probes provide a rough estimate ofbackground based on coarsely modeled sequence, without taking intoaccount any spatial variation. Use of BGP over MM probes requiresroughly 50% less space on the array; thus use of BGP probes allowsroughly a 2 fold increase in PM content.

Background estimates from BGP and MM probes may not be accurate forspecific PM probes leading to biases in the background estimatesresults. We therefore look at the typical distribution of error inmodeling background, as well as the downstream results on the ability todetect differential change. These results are based on a research arraydesign using a pre-optimized version of the whole transcript assay.

We examine the performance of these two methods with a controlledexperiment on a research chip using spiked-in transcripts. Thisexperiment consists of 9 full-length clones which are spiked in atseveral different levels of relative abundance, with three replicatesper level of abundance. The levels are arranged in a standard latinsquare design so that each experimental pool has some transcripts ateach level of abundance, and every transcript occurs at all levels ofabundance. In this experiment, the standard HeLa background is used as acomplex background, and is not varied. This experiment was performed onan initial testing array design (“WTA-Sensor”) with a pre-optimizedversion of the whole transcript assay. protocol. This array designcontained probes designed to complement every other base in thefull-length transcripts. For analysis purposes, the transcripts weredivided into twenty sub-regions each (simulating exons), and four probeswere chosen to represent each sub-region, leading to twenty probe setsper transcript, or 180 probe sets analyzed in the data set used here.

In FIG. 5, the distribution of multiplicative error for three differentbackground models is displayed for a collection of background probes. Bydesign, these should all be without specific hybridization, and hence anideal background model would exactly predict the intensity of theseprobes. What is displayed in FIG. 5 is a boxplot of the residualslog(I)−log(B) for each model. The first model is the log-scale errorinvolved in approximating the background for probes of differing GCcontent by a single, universal background quantity, in this case, themedian intensity. The second model sets the background estimate for agiven probe to be the median of the corresponding bin of probes withsimilar GC-content, i.e log(I)−log(B(GC)). The last model sets thebackground estimate to be the corresponding MM probe. Both MM and GC-binmedians closely follow the behavior of probes across a wide range ofGC-content, although GC-bin medians are slightly biased. Using a singleglobal estimate of background has an unacceptably large distribution oferrors in background prediction, however, GC-bin medians aresufficiently low error to be considered as replacements for specific MMprobes.

FIG. 6 shows the mean-squared log-scale error resulting from each of thethree considered models. Use of a single, global background estimateleads to large mean-square-error, and so is not acceptable forbackground estimation. Despite the very similar inter-quartile rangesfor both the gc-bin background model residuals and the mm modelresiduals, the mean-square error is larger for gc-bin probes. However,the known tendency of mismatched probes to hybridize to the true targetsuggests that gc-bin background estimation may be sufficient forbackground estimation in the context of detecting differential changesin signal.

Using the same data set, performance was evaluated for the variousbackground estimation methods. In this case, performance is measured bythe ability to detect differential change between paired experiments. Inthe figure below we show receiver operator curves (ROCs) based onthresholding a t-like statistic between the different spikeconcentrations. Use of BGP over MM probes does not adversely affect theability to detect differential change. It was observed however that someprobesets exhibited more positive bias with BGP, indicating thatcoarsely modeled background with BGP does not completely remove biasfrom the data to the extent of the PM-MM method. It is important to notethat these results are based on an assay that generates DNA target, notRNA target. Thus the observations for the whole transcript assay withregard to BGP may not hold for the 3′ IVT assay.

Genomic and Antigenomic Background Probes

Genomic Background Probes: Background probes were selected from aresearch prototype human exon array design based on NCBI build 31.Mismatch probes were designed against Genscan Suboptimal exonpredictions that were not supported by any other annotation source. Whenpossible, 1000 mismatch probes were selected for each bin of GC contentfrom 0 to 25.

Antigenomic Background Probes: The Antigenomic background probesequences are derived from 15-mer root sequences that are not found inthe human (NCBI build 34), mouse (NCBI build 32), or rat (HGSC build3.1) genomes. These root sequences were extended by 5 bases on each end,where the added bases were drawn randomly from a GC distributionidentical to the base 15-mer. There are probes for GC counts rangingfrom 2 to 24 per probe sequence. When possible, 1000 mismatch probeswere selected for each bin of GC content.

Exon Probeset Annotations and Transcript Cluster Groupings

Procedures for grouping and annotating probesets are also disclosed.Appropriate grouping of probesets into transcript clusters andsubsequent filtering of probesets within the transcript cluster plays acritical role in generating gene-level signal estimates. The annotationsand groupings provided are intended to be a baseline of information foreach exon array. How the information is applied will depend on thespecific goals for each experiment. Definitions for terminology andarray design are disclosed in the GeneChip Human Exon 1.0 ST ArrayDesign Technote which is incorporated herein in its entirety for allpurposes.

The GeneChip Human Exon 1.0 ST Array and other exon arrays are moreexploratory than predecessor expression arrays such as the HG-U133 2.0Plus. Previous array designs used a cDNA assembly consensus or exemplarapproach where one or a few probeset were associated with a specificgene. For the exon arrays the design was exon focused rather than geneor transcript focused. As a result there are no intrinsic transcript orgene entities in the exon array designs. Instead there are onlyprobesets associated with exons or contiguous parts of exons. Groupingsof exon probesets into transcripts and genes is now a dynamicpost-design process.

As new genome assemblies and annotations are released there is anopportunity to generated improved transcript cluster groupings. Themeta-probeset lists available in the support files for the exon arraysare one such example of these groupings. Meta-probeset lists fromdifferent versions of the genome are likely to have different groupingsof exon probesets into transcript cluster. Furthermore these differentgroupings may in turn have subtle (or even substantial) impacts on genelevel analysis.

The basic approach used to generate transcript cluster groupings for aparticular genome is to (1) construct gene annotations on the genomeusing a variety of annotation sources merged using a set of rules, (2)map exon probesets to gene annotations using the genome, and (3) use thetarget genome to group the probesets that mapped to the geneannotations. A fall out from this is that gene annotations are easilycreated from the transcript annotations used to group the probesets. Thegene grouping process and probeset annotation process are described inmore detail below.

The result of the entire process is a collection of “Design AnnotationFiles”. These design annotation files are available from the appropriateexon array support page on the Affymetrix web site. Specifically, foreach genome version a set of master GFF files are generated. These filesare used to populate parts of the NetAffx content and are used togenerate the other files including: design level probeset annotation CSVfile, meta-probeset list files and IGB binary annotation (BGN) files

Defining Gene Annotations

The following five steps outline one method for grouping probesets. (1)Cluster transcript annotations on the same strand of the target genomeusing a set of rules involving exon overlap and splice site sharing. (2)Label each transcript cluster as a different gene. (3) Join exons ofclustered transcript annotations on the target genome to determine genestructure. (4) Map each probeset to a single gene, based on whether itfalls within the gene's annotated exon boundaries on the target genome.(5) Group probesets together that map to the same gene as a transcriptcluster.

There are situations where simply clustering transcripts based onexon-exon overlap can produce gene annotations that actually representmultiple genes. Sometimes this is probably due to two genes actuallysharing some transcribed region on the genome (i.e. 3′ UTR of one geneoverlapping the 5′ UTR of a downstream gene on the same strand). Inother cases this is due to erroneous cDNA sequences, alignmentalgorithms, or gene predictions. There was an effort to avoid suchover-clustering of transcripts, even at the expense of fragmenting somegene annotations. This decision was motivated by the fact that falselyjoined exons (from more than one gene) would generate results which atfirst pass would appear to be alternative splicing.

Gene annotations were constructed from clusters of transcriptannotations. Because of the wide variety of transcript annotationssources, it was desirable to find a way to qualitatively categorize thedifferent types of sources. For example, RefSeq sequences are manuallycurated and generally trusted as accurate isoforms of genes. While ESTsand partial mRNAs are not manually curated, they cover greater portionsof the target genomes and are supported by the fact that they arederived from cloned sequence. And although ab-initio gene predictionsare not inherently supported by clone sequence, they may be reasonablepredictors of transcript structures. Acknowledging that differentannotation sources have different levels of confidence, an iterativeapproach was implemented to ensure that high confident annotations werenot merged by less confident annotations.

For the purposes of establishing a hierarchy of gene confidence levels,the sources of input transcript annotations were partitioned into threetypes. From highest to lowest confidence, the types were labeled core,extended, and full. Broadly defined, the core type consisted of (BLAT)alignments of mRNA with annotated full-length CDS regions, the extendedtype consisted of cDNA alignments and annotations based on cDNAalignments, and the full type consisted of sets of ab-initio genepredictions.

For the GeneChip Human Exon 1.0 ST Array, the annotation sources foreach level are: Core Gene Annotation sources, RefSeq alignments, Genbankalignments of ‘complete CDS’ transcripts, Extended Gene Annotationsources, cDNA alignments, Ensembl annotations (Hubbard, T. et al.),Mapped syntenic mRNA from rat and mouse, microRNA annotations, Mitomapannotations, Vegagene (The HAVANA group, Hillier et al., Heilig et al.),VegaPseudogene (The HAVANA group, Hillier et al., Heilig et al.), FullGene Annotations, Geneid (Grup de Recerca en Informatica Biomedica),Genscan (Burge, C. et al.), GENSCAN Suboptimal (Burge, C. et al.),Exoniphy (Siepel et al.), RNAgene (Sean Eddy Lab), and SgpGene (Grup deRecerca en Informatica Biomedica), TWINSCAN (Korf, I. et al.).

The core type was so named because the annotations in this type wereintended to be the foundation from which we built our gene annotations.The extended type derived its name from the sense that these annotationswould extend the boundaries of the core genes. The idea behind the nameof the full type was that it would signify all possible content.

The transcript annotation clustering procedure utilized the confidencerankings of the transcripts in such a way that non-overlappingtranscripts from a higher rank could not be joined together by atranscript annotation from a lower rank. For example, annotations of ESTalignments would not be used as evidence to join two separate RefSeqalignment annotations—however, a third RefSeq alignment annotationcould. Lower ranking annotations could only be added to the transcriptclusters established at a higher ranking, or they would establish newclusters without any higher ranking content.

In one embodiment the algorithm for clustering the transcriptannotations is as follows:

-   1) Core annotations are merged first.    -   a. Spliced core annotations are merged first.        -   i. If two core annotations share a splice site, they belong            to the same cluster.        -   ii. If two core annotations overlap with two different exons            in each transcript, they belong to the same cluster.        -   iii. If a core annotation lies within the boundaries of            exactly one other core exon, the two belong to the same            cluster.    -   b. Single exon core annotations are added next.        -   i. If a single exon core annotation overlaps an exon of            exactly one spliced core cluster, it is added to that            cluster. (However, the single exon is not used to determine            exon overlap for subsequent single exon core            annotations—i.e. prevent “chaining” single exons).        -   ii. If a single exon core annotation lies within the            boundaries of exactly one other core exon, the two belong to            the same cluster.        -   iii. Single exon core annotations that don't overlap any            spliced core annotations are merged with each other based on            overlap.        -   iv. If a cluster of single exon core annotations formed in            step (iii) overlaps with exactly one other spliced core            cluster (with the single exons from step (i) now included),            it is added to the spliced cluster (i.e. now allow chaining            if it is completely unambiguous)-   2) Extended annotations are added next, in a similar manner as the    single exon core annotations.    -   c. If an extended annotation exon overlaps the exon of exactly        one core cluster or shares a splice site with exactly one core        cluster, it is added to that cluster. (However, the added        extended annotation is not used to determine exon overlap for        subsequent extended annotations.)    -   d. If an extended annotation falls entirely within the bounds of        exactly one core exon, the extended annotation is added to that        core cluster.    -   e. If a spliced extended annotation overlaps the exons of more        than one core cluster and does not share a splice site with        exactly one of them, the extended annotation is broken up into        its underlying exons. These underlying exons are then treated as        separate single exon transcripts for the purposes of merging.    -   f. Extended annotations that don't overlap any core clusters are        merged with each other.    -   g. If a cluster of extended annotations formed in step (d)        overlaps exactly one other core cluster (with the extended        annotations from step (a) now included), it is added to that        core cluster.-   3) Full annotations are added last, in exactly the same manner as    the extended annotations, except that extended clusters are treated    the same way as core clusters.    FIGS. 7-10 illustrate hypothetical examples of how some of the rules    for clustering transcript annotations may be applied.

The resulting transcript clusters define a set of exon boundaries, whichin turn determine the gene definitions used for probeset clustering.Each gene annotation is constructed from transcript annotations from oneor more confidence levels. Some parts of a gene annotation may derivefrom high confidence core annotations, while other parts derive from thelower confidence extended or full annotations. Therefore, differentparts of the gene annotation can be labeled according to the highestconfidence level of transcript annotation that supports that part. Thislabeling of the resulting gene annotations according to confidence levelwas used to further annotate the probesets that mapped to the gene.

Probeset Gene Mapping and Probeset Grouping

Probesets were grouped together if they mapped to the same geneannotation. Generally, a probeset was said to map to a gene annotationif it fell inside the exon of that gene. However there were cases thatwere not so straight-forward, and so a set of specific rules wasestablished. (1) If a Probe Selection Region (PSR) overlaps exactly onegene, then the associated probeset is mapped to that gene (even if theprobeset does not overlap any genes). (2) If a probeset falls entirelywithin the exon bounds of exactly one gene, then it is mapped to thatgene. (3) If a probeset falls entirely within the exons of multiplegenes, the probeset is mapped to its own unique transcript, unless theprobeset falls within the exons of multiple genes, but within the coreregion of only one gene, then it is mapped to the gene with the coreregion (4) If a PSR or probeset does not overlap any genes, the probesetis mapped to its own unique transcript.

There were some cases where probesets were selected from Probe SelectionRegions (PSRs) that were based on content that was removed after thedesign. Rule 1 attempted to associate the probesets from these PSRs withthe genes that they were proximal to.

Cases where probesets fell within exons from multiple genes were coveredby Rule 3, however there were special cases involving core genes. If aprobeset fell within the exons of multiple genes, but only one of themwas a core gene, then the probeset was mapped to the core gene. Themotivation behind this rule was that since core genes were regarded withthe highest confidence, a probeset that fell within core exon should bemapped to that gene, despite the lower confidence content.

There were three types of annotation that were applied after theprobesets had been grouped: (1) a confidence ranking, (2) whether theprobeset was ‘bounded’ to a transcript cluster (less confidentassociation) or actually contained (high confident association) and (3)whether the probeset was within a high quality CDS region. While allprobesets were given a confidence ranking, the ‘bounded’ and ‘CDS’ flagswere applied to only some probesets.

The annotations used to generate the gene bound definitions were furtherused to rank the probesets with respect to their confidence level. Asmentioned above, the gene annotations used to group the probesets hadcore, extended, or full regions according to the confidence level of thetranscript annotations used to compose the gene. A probeset with aposition in the gene that was within a core region of the gene was givena ‘core’ level ranking. Likewise for probesets that fell within theextended and full regions of their associated genes. In the cases wherea probeset crossed level boundaries, the probeset was labeled with thelower confidence rating. If a probeset did not overlap any genes, it waslabeled as ‘free’. See FIG. 11.

There were other instances where a probeset fit inside the exon of morethan one gene annotation, and no determination could be made as to whichgene annotation it belonged to. In these cases, these probeset waslabeled ‘ambiguous’ and placed singly in its own gene annotation. Theexception to this rule was if a probeset mapped to the exons of multiplegenes, but fell within the core region of only one gene. Then theprobeset was given the core level annotation instead of ‘ambiguous’.

Many of the genes generated by the transcript clustering procedure wereactually single exon genes defined by either a solo GENSCAN Suboptimalexon annotation or a single exon EST alignment. These annotationsoccurred at relatively high frequency throughout the greater transcribedregions of the target genome. Many of these annotations occurred withinthe introns of larger, spliced, higher confidence genes. It was decidedto include these probesets with the transcript clusters they fall in,but to provide annotations to the effect that the grouping confidence islower. For example, this allows researchers to easily include or excludethese single exon probesets for alternative transcript analysis with thesame locus.

Therefore, if a probeset mapped to a single exon extended or full genethat lay within the intron of exactly one other gene, then the probesetwould be “regrouped” as part of the larger, spliced gene and given theadditional label, ‘bounded’. The intuition behind the word ‘bounded’ isthat the probeset is bounded by the larger spliced gene. See FIG. 12.

Probesets that fell completely within the CDS region of at least onetranscript from one of the following four sources were flagged as ‘cds’probesets: RefSeq alignments, Genbank alignments ‘complete CDS’transcripts, Ensembl annotations, and VegaGene annotations. This allowsresearch to quickly focus down to a set of probesets against likelyprotein coding regions.

The probesets were further characterized by listing the quantity andsource of annotations that could be interrogated by the probeset. Thesame transcript annotations used as building blocks for gene definitionswere also used to describe the evidence for the probesets (although, ingeneral any annotation set can be used here). As with the other labelingprocedures, the entire probeset had to be bounded by the exons of theannotation in order for it to be listed as evidence. This technique forlisting evidence gave rise to some non-intuitive situations where aprobeset could have been given a confidence level of ‘core’, ‘extended’or ‘full’, but not have any evidence. This is due to the fact that thetranscript annotations were merged into larger gene annotations whendetermining the confidence levels, but were not merged when determiningevidence for the probesets. Therefore, a probeset could be mapped to anexon that was derived from a composition of several smaller annotations,yet the probeset did not fall within the bounds of any one of thesesmaller annotations. See FIG. 13.

The mapping of exon array probesets to specific annotations is containedwithin the GFF file as comments. These GFF comments can be easily parsedshould this level of information be needed.

These probeset groupings and annotations are meant to provide a startingpoint for researchers and are not intended as the final word on howprobesets should be clustered and analyzed. While the notion oftranscript clusters loosely equals a gene, it should be noted that wemake no attempt to merge non-overlapping transcript clusters for thesame gene (based on paired EST reads for instance).

Gene Signal Estimates from Exon Arrays

With exon arrays like the GeneChip® Human Exon 1.0 ST Array, researcherscan examine the transcriptional profile of an entire gene. Being able togather data for each individual exon enables the investigation ofphenomena such as alternative splicing, alternative promoter usage, andalternative termination while also providing more probe level data todetermine the overall rate of expression a particular locus. Thecalculation of a gene level signal estimate is of particular interestfor many analysis as it provides a familiar starting point for theanalysis of expression data. Generating gene level signal estimates fromexon array data is the focus of this whitepaper. Throughout the rest ofthis article, the phrase “transcriptional locus”, or “locus”, willdescribe a genome location that is transcribed rather than using theterm “gene” which often implies protein-coding transcripts.

The GeneChip® Human Exon 1.0 ST Array is inclusive by design. Inaddition to containing probe sets that interrogate exons of RefSeq(Pruitt and Maglott, Nucleic Acids Res 29: 137-140, 2001) genes, mRNAsand ESTs from GenBank (Benson et al. GenBank. Nucleic Acids Res 27:12-171999), it also contains probe sets for exons predicted by ab-initiogene finders such as GENSCAN (Burge and Karlin J Mol Biol 268: 78-94,1997), TWINSCAN (Korf et al. Bioinformatics 17 Suppl 1: S140-148, 2001),geneid (Parra et al. Genome Res 10: 511-515, 2000), and even GENSCANSuboptimal exon predictions. Having so much content on the array allowsthe detection of novel alternative events, but presents a new challenge.This is because including alternative exons can negatively affect theestimate for the entire transcriptional unit. Additionally, much of thecontent is exploratory in nature, such as the GENSCAN suboptimalpredictions, and by nature will have a lower probability of beingactually transcribed compared with other well-annotated sequences.Having such a rich set of possible probe sets for each locus providesboth the benefit of more data and the challenge of finding the bestprobe sets to use for gene-level signal estimation in a givenexperimental study.

The first step in estimates signal for a particular locus is todetermine the boundaries for each locus, or gene, and which probe setsare contained within the boundaries. It is important to avoid mistakenlyjoining two loci, or genes, that are transcribed separately or splittinga single locus in two or more loci. Affymetrix has used genomicannotations such as RefSeq, cDNAs, and gene predictions to creatediscrete locus definitions (Wheeler, Gene Bounds: Probe Set Grouping,Annotation, and Evidence. In Affymetrix White Papers.2005). Briefly,annotations that overlap exons on the same genomic strand wereiteratively merged into a larger locus annotation. The probe sets on theGeneChip Exon Array were then assigned to these locus annotations, thuspartitioning them into distinct transcriptional units, corresponding inmany cases to genes.

Different, sometimes nested, loci definitions were created for differentlevels of annotation confidence so that research can clearly identifymore exploratory content from well annotated content. High confidenceannotations such as RefSeq transcripts and full-length mRNAs make up the“core” loci. The core annotations in addition with all other cDNA basedannotations make up the “extended” loci. The “full” loci are theextended loci plus the ab-initio gene predictions. There are many singleexon annotations that do not overlap the exons of other multi-exonannotations, but are contained in the intron of a larger multi-exonannotation. In some cases, such as a suspected alternative exon, it maybe desirable to associate these single exon annotations with the locusthey are nested within. This is referred to as being “bounded” by anannotation.

Users may select what confidence level of loci definitions to use andalso exclude probe sets associated with a particular annotation. Forexample, the exploratory GENSCAN Suboptimal probe sets will not beappropriate for studies examining transcription changes inwell-annotated genes.

To explore the effect of uncertainty of which probe sets to use forestimating the signal of a particular locus using PLIER (Hubbell, PlierWhite Paper. In Affymetix White Paper., 2005) signal estimates, thefollowing experiments wereerformed:

-   -   Addition of decoy probe sets likely to have low signal to        simulate the inclusion of speculative content that is rarely, if        ever, transcribed.    -   Addition of decoy probe sets likely to have high signals that do        not correlate with the original locus to simulate the effect of        mistakenly adding probe sets from a transcriptionally distinct        locus.    -   Effect of using probe sets from more speculative annotations.    -   Effect of adding probe sets likely to have low signals on the        signal estimates of spike-in data sets.

Addition of Decoy Probe Sets:

Rather than generating simulated data we will sample from the GENSCANsuboptimal only probe sets as a proxy for low signal probe sets (meanPLIER estimated signal across all GeneChip arrays: 76+/−574) and thefull length mRNA probe sets for high signal probe sets (mean PLIERestimated signal across all GeneChip arrays: 339+/−2,654). To facilitatethe study of adding low- and high-signal decoy probe sets, we identifieda set of probe sets that appeared to be constitutive for 1674 loci, with4-11 constitutive probe sets in each locus, on chromosomes 1, 10, 11,12, and 22. These constitutive probe sets were used as a gold standardto generate PLIER signal estimates for those loci. For this experiment,we used a set 48 GeneChip® Human Exon 1.0 ST Arrays hybridized to 16different tissues with 3 technical (assay) replicates for each tissue.Tissue replicates were quantile normalized and then the differenttissues were median normalized to each other. For both the low- andhigh-signal decoy probe sets, we then generated pair wise correlationsfor each PLIER estimate of the test case vs the PLIER estimate for thegold standard set (FIGS. 14 and 15).

PLIER appeared relatively robust to the addition of limited amounts oflow-signal decoy probe sets (FIG. 14), but was more affected by theaddition of higher-signal decoy probe sets (FIG. 15). It should be notedthat the addition of low-signal decoy probe sets will also reduce theabsolute values of the PLIER estimates. As feature response of probes isconstrained to have an average of 1, the low-signal probe sets make theother probe sets look like they are more sensitive to the presence ofthe target in comparison and thus reduce the signal estimate (seeHubbell 2005 for discussion of constraints)

Using More Speculative Annotations:

In addition to adding artificial decoy sets to the gold standard setdescribed above, we also observed the effect of adding probe setsresulting from annotations that are more exploratory in nature for aparticular locus. As probe sets with more speculative annotations areincluded in the loci expression estimates, the correlation with abioinformatically enriched constitutive set becomes worse. Theannotation file provided with Human Exon array details the annotationsassociated with each probe set. Using just the probe sets resulting fromfull-length mRNA and RefSeq exons results in loci expression estimatesthat are very highly correlated.

Expression estimates from extended annotations, those resulting from anycDNA, are also highly correlated with the gold standard set. Expressionestimates from full annotations, those resulting from a gene predictionprogram, were also highly correlated, but not as much as the cDNA basedannotations.

Including probe sets that were bounded by a gene (located in an intron,but not overlapping an exon) resulted in the lowest rate of correlation.The effect on a particular gene can be seen in the form of both lowerPLIER estimates and lower correlation with the constitutive set across48 samples. As more exploratory probe sets are added to thetranscriptional locus estimates, the overall correlation is still high,but the range of PLIER estimates becomes compressed.

Effect of Decoy Sets on Differentiating Known mRNA Concentrations:

We also investigated the effect of adding low-signal probe sets on theability to differentiate between spiked-in concentrations of mRNAs. Forthis experiment, a latin square design was used with 4 groups ofspike-in RNAs at 4 different concentrations of 0, 0.18, 0.36, and 0.71pM with 3 replicates each. A receiver operator curve (ROC) wascalculated with the ability to differentiate between spike-ins that arepresent in different concentrations and those present at the sameconcentrations (FIG. 16). The addition of low-signal decoy probe setsdecreases the ability of the PLIER expression estimates to differentiatebetween spike-in concentration pools (FIG. 17).

Strategies for Coping with Exploratory Content:

One natural solution to determine which probe sets to use forsummarizing a particular locus is to only use probe sets containedwithin high quality annotations such as RefSeq transcripts. Leveragingyears of work cloning and sequencing genes combined with expert manualannotation assures that the probe sets chosen belong to that locus. Thisapproach does not address the issue of including probe sets toalternative events, but if most exons are constitutive for a given gene,they should dominate the PLIER expression estimate rather than thealternative probe sets. Robust methods such as PLIER and RMA (Irizarryet al. Nucleic Acids Res 31: e15, 2003) should be minimally affected bya limited amount of alternative splicing at a particular locus.

Selecting Features to Use in a Given Data Set:

Another strategy is to use feature sets from more exploratoryannotations and iteratively discard those that appear to be performingpoorly. This approach takes advantage of PLIER's ability to identifysome of the signal at a particular locus and iteratively excludefeatures that are not correlated with that signal. As aproof-of-principle study, we implemented a procedure, IterPLIER, to useall features to generate an initial PLIER estimate and found the 22features best correlated with that estimate, re-estimated the expressionestimate with those 22 features using PLIER and then selected 11features that were best correlated with the new expression estimates togenerate a final PLIER signal estimate.

The correlation with the signal estimates without decoy probe setscompared to the result of IterPLIER on the low-intensity decoy set (FIG.18) and high-intensity decoy set (FIG. 19) are improved to thoseresulting from using the initial PLIER estimates. The results forIterPLIER on the Akt3 locus illustrate that signal estimates consistentwith the original bioinformatically enriched constitutive gold standardset were obtained for all cases but the full-bounded set. These resultscompare favorably to the raw PLIER expression estimates. Thisillustrates the improved correlation for all but the full-bounded probeset groups. This method also improves the ability of PLIER expressionestimates to differentiate between spike-in concentration pools in thelatin square experiments. See FIG. 20.

Determining the best way to estimate locus intensity, and hence theunderlying relative target response, is still an area of activeresearch. Here we have shown that PLIER is more robust with respect toadditional low-signal probe sets than to additional conflictinghigh-intensity probe sets. These results suggest that loci which may betranscribed separately should be joined with caution, while lowconfidence annotations may be included in an analysis with less risk. Asthe new Human Exon array has more data for each locus we now have anopportunity to improve our estimates, but also face more choices aboutthe type of estimates that are most appropriate for the study at hand.

Quality Assessment of Exon Arrays

Quality assessment can play an important role at many stages in theprocess of generating gene expression data for biological samples, fromarray design to sample processing to signal estimation and analysis. Inthis report we describe some quality assessment procedures that arebased on CEL intensity data from the GeneChip® Human Exon 1.0 ST Array.The methods detailed here are described in Chapter 3 of the Bioconductormonograph [Robert Gentleman, Vince Carey, Wolfgang Huber, RafaelIrizarry, Sandrine Dudoit, Bioinformatics and Computational BiologySolutions using R and Bioconductor, 2005.]. The quality assessmentprocedures we consider entail computing some summary statistics for eacharray in a comparable set and comparing the level of the summarystatistics across the arrays. We therefore assume that the user has aset of comparable arrays of the sort that would normally be analyzedtogether to address substantive biological questions. The qualityassessment part of the procedure entails identifying outliers within theset. These could be flagged and excluded from the substantive biologicalanalysis, or the downstream analysis could be adjusted to account of theoutlier arrays, by weighting for example.

The procedures we describe can identify outliers, but it does notprovide hard and fast rules (specific thresholds) as to which arrays toflag as outliers. These need to be developed in the context ofparticular applications with specific types of samples, with knowledgeof the costs involved in repeating experiments and the costs of drawingwrong conclusions. The procedures also do not provide guidelines toassess the overall quality of a set of arrays. These can be arrived atby collecting data for many sets of arrays and using them to set targetlevels and variance limits for the various quality assessment measureson a user-by-user basis. Lastly, while the methods presented here can beused to identify outlier arrays, they will not indicate why the array isan outlier (ie input RNA quality problem, target prep problem,hybridization problem, etc . . . )

The following are included below: simple summaries of the intensitydistributions and log₂ ratios of cel intensities, summaries of featureset model fit derivates—signal estimates and residuals, and the specificquality metrics reported in the Exon Array Computational Tool (ExACT)Quality Report.

Feature or cel level summaries. Note that the feature level graphspresented below are not directly obtainable from the Exon ArrayComputational Tool (ExACT). The “exact-probe-intensity-extraction”command line program can be used to extract probe level intensities foreither the whole array or a subset of the features on the array.Normalized and un-normalized CEL files can be used with this tool. Theextracted probe intensities can then be evaluated with a statisticalapplication such as R, S-Plus, or MatLab to generate the summary graphs.

Single Array Cel Intensity Distribution Summaries

The simplest and most obvious thing to do with a set of cel files thathave been collected for the purpose of addressing some questions ofbiological interest is to examine the distribution of cel intensitiescorresponding to each array in the set. As the array distribution ofintensities is highly skewed, the log₂ transformation helps toapproximately normalize the distributions. Log base 2 is typically usedas it facilitates the conversion back to the original scale: adifference of 1 on the log base 2 scale corresponds to a fold change of2 on the original cel intensity scale. Array log₂ cel intensitydistributions can be summarized by a three point summary: 25^(th)percentile (Q1), 50^(th) percentile or median (Q2) and 75^(th)percentile (Q3). To identify outlier arrays within a set, it helps tocompare the distribution by means of boxplots. In the boxplot display, abox is formed with sides at the 25^(th) and 75^(th) percentiles of thedistribution. A line is also drawn within the box at the level of themedian. Whiskers are also drawn extending beyond each end of the boxwith points beyond the whiskers typically indicating outliers. As thepurpose of the boxplots in this application is to compare the locationand spread of the main body of the distribution of log₂ cel intensities,we typically suppress the plotting of symbols for points beyond thewhiskers. See, Tukey, John W., Exploratory Data Analysis,Addison-Wesley, Reading, Mass., (1977) for guidelines regarding outlierdetection.

FIG. 21 shows an example set of boxplots of log₂ cel intensities for aset of 20 arrays. The 20 arrays were hybridized with preparations madefrom RNA extracted from colon tissue samples. We observe a range of celintensity distributions across the 20 arrays. A normalizationtransformation is typically applied to make the distribution of celintensities more comparable. Based on cel intensity distributions aloneit is difficult to assess the impact that observed differences will haveon downstream analysis. It is recommended that cel intensitydistributions be observed and recorded for future reference.

Log₂ intensity distributions can also be summarized by means of ahistogram display. Whereas histograms provide more detail, but makes thecomparison of several distributions more difficult. Histograms enablethe detection of a secondary mode in the distribution corresponding to abright spot on the array. Unless the array images are examined toidentify locally bright or dim regions, it may be warranted to examinethe histograms of the log₂ cel intensity distributions for any evidenceof bright spots. FIG. 22 shows the histogram of log₂ cel intensities forthe first 8 arrays in the set of 20 displayed in FIG. 21. As there is noevidence of bi-modality, these displays provide no additionalinformation to what is captured by the boxplots in FIG. 21.

Multiple Array Feature Level Summaries

A significant characteristic of CEL intensity data is the largemagnitude of feature specific effects. The typical range of log₂ CELintensities for arrays in FIG. 21 is (lo, hi) (e.g. (2, 16)). The rangeof variability for any given feature, or CEL, across the set of arrayswill be much lower—bright features tend to be bright across all arrays,and dim features tend to be dim across all arrays. The model used byPLIER (see PLIER Technote) to compute signal estimates, as well as themodel used by RMA to compute feature set summaries, account for thesefeature effects. The MA-plot [Dudoit S, Yang Y H, Callow M J and Speed TP. Statistical methods for identifying differentially expressed genes inreplicated cDNA microarray experiements. Statistica Sinica, 12:111-139,2002.] provides a visual display that corrects for feature-specificeffects and facilitates making comparative assessments of cell intensitydistributions that are sensitive to the technical variability which wewant to detect. Suppose that x and y are two vectors of log₂ intensitiescorresponding to two arrays. The MA-plot is a scatter plot ofM_(i)=(x_(i)−y_(i)) on the vertical axis vs. A_(i)=(x_(i)+y_(i))/2 onthe horizontal axis. The M values are relative log₂ intensities, or log₂ratios of intensities. The median M value gives a sense of relativeshift in log₂ intensity between two arrays, and the inter-quartile range(Q3−Q1) of M values gives a sense of reproducibility of cell intensityvalues between arrays. The M values also contain some real biologicalvariability, but for practical purposes that variability in M valuesprovides a useful indicator to detect extra technical variability. TheMA-plot as described above provides pair wise comparisons between anypair of arrays. When assessing quality in a set of arrays, pair wisecomparisons produce much redundancy. To avoid this redundancy, eacharray can instead be compared to a synthetic array created by takingfeature-wise medians across the set. For example application andinterpretations of MA-plots see Dudoit et. al.

Boxplots of M values provide a summary of MA-plots which ignores therelationship between relative log₂ intensity and mean log₂ intensity butfacilitates the comparative assessments of the distribution M valuesamong many chips. FIG. 23 shows boxplots of relative log₂ intensityvalues for the 20 arrays whose log₂ intensity values are displayed inFIG. 21. The two figures, FIG. 21 and FIG. 23, deliver the same message,with latter being more sensitive to differences between arrays. At thispoint, a word on normalization is in order. Ideally, the distributionsof log₂ cell intensities for a set of arrays to be analyzed togetherwill be comparable—the three-point summary values (Q1, Q2, Q3) should bearound the same level for all arrays, and the MA-plots should be flat,narrow and centered at zero. This ideal will rarely be achieved inpractice due to the biological and technical variability that is aninherent characteristic of the cell intensity data being produced by avery elaborate process with multiple sources of variability.Normalization is the term used to refer to procedures that adjust cellintensity data in order to remove the technical variability withoutmasking the biological variability. Some normalization methods willremove much of the variability that is detected by the qualityassessment displays describe above. Assessing variability in theun-normalized data values is still important as the amount of adjustmentthat is required to normalize the cell intensities of an array is a goodindicator of quality—the less adjustment required, the better thequality.

Summaries of feature set model derivatives. As mentioned above,feature-specific effects account for a large fraction of the variabilityin intensity between features. These feature effects are a function ofsequence dependent probe affinities as well as the composition of thehybridization mixture hybed to the array. Models have been developed toaccount for feature effects by estimating these effects in a multi-arraycontext. The models are typically fitted to a biologically meaningfulset of arrays—RNA extracted from liver samples; samples characterized asnormal or cancerous, for example. In this section we discuss qualityassessments based on derivatives, or by-products, of the PLIER modelfits. The models are fitted feature set by feature set and producefeature effect estimates, signal estimates and residuals. Featureeffects are a nuisance parameter which yield improved signal estimatesbut can otherwise be ignored. Signal estimates are the expressionindicators that will be used in the biological applications. They canalso be used for quality assessment, as can the model residuals. Thenext two sections describe quality assessment measures that are based onthese model-derived quantities. The analysis presented below can beperformed in a variety of statistical analysis packages using the ExACTprobeset summary output files as well as the option PLIER residualfiles.

Relative log₂ signal summaries. After fitting the PLIER model to a setof arrays, we have signal estimates for each feature set and each array.We can assess data quality in terms of reproducibility of signal valuesacross arrays as measured by log₂ ratios of signal estimates or relativelog₂ signal. In a set of arrays hybridized to samples coming fromdistinct biological samples, some signal variability is expected due toreal biological variability. Using a measure of variability based on themost reproducible relative log₂ signal values—an inter-quartile rangewhich in effect uses the 50% most reproducible signal values to assessreproducibility—produces a measure which can detect extra technicalvariability above and beyond biological variability. We can assess thereproducibility of signal estimates between two chips by means of anMA-plot, in which relative log signal is plotted on the vertical axisand average signal on the horizontal axis. As pair wise comparisonsproduce much redundancy, signal values for each array can instead becompared to signal values from a synthetic array constructed bycomputing for each feature set the median signal over all the arrays inthe set. Superimposing the MA-plots with smooth curves capturing thequartiles (median, lower and upper quartiles) in relative log₂ signalvalues as a function of average log₂ signal is a helpful visual aid todetect subtle differences among MA-plots.

Boxplots of relative log₂ signal values provide a succinct summary ofMA-plots which ignores the relationship between relative log₂ signal andmean log₂ signal but facilitates the comparative assessments of thedistribution relative log₂ signal values among many chips. FIG. 24 showsboxplots of relative log plier expression values for the set of arraysanalyzed in FIG. 21. We see that with the except on of array number 8,at the plier signal level the arrays are comparable in expressionlevel—the center of the boxplots close to zero—and reproducibility—thesize of the boxplots are comparable. Based on these assessments, onewould be advised to exclude the data from array 8 from downstreamanalyses.

In the assessment of reproducibility discussed so far we have implicitlyassumed that the summaries would be computed based on all feature sets.We can generalize the above assessments by summarizing over a variouscollections of feature sets. One may wish, for example, to look atsignal reproducibility in the normalizing exon feature sets, or thenormalizing intron feature sets, or a random sample of feature sets.

Summaries of residuals from fits. A by-product of the PLIER fit is theresiduals from the model fit. These can be summarized in various ways toproduce quality indicators. The median of the absolute value ofresiduals, pooled over all features sets, provides a simple array-levelsummary that will detect some differences in technical variabilitybetween arrays. FIG. 25 displays a bar graph, with bar heightproportional to mad(plier residual) for the 20 chips we have been usingfor illustration. Array 8 stands out as the chip with the largestmad(plier residual). Mad(plier residual) is highly correlated with theinterquartile range of relative log₂ signal.

Alternatively, the residuals can be combined into standard errors ofsignal estimates, and these can be summarized over all features sets onthe array. The standard error can be normalized, by the median standarderror for each feature set, for example, to provide a quality indicatorthat is slightly more sensitive to differences between arrays (seederivation of NUSE values in the Bioconductor monograph). As in the caseof the relative log₂ signal values summaries, the various residualsummaries can be confined to various collections of features or featuresets.

ExACT Quality Report. The probe summarization step(exact-probe-summarize command) will generate an optional report file.This text file contains various summary metrics from the CEL file data.It is important to note that the report file results are dependent onthe particular analysis parameters used; thus different results can begenerated from the same set of CEL files depending on which intensitymethod, quantification method(s), and probesets are used.

The quality report file includes meta information about the particularanalysis run. This information is encoded as lines starting with thepound symbol (‘#’). For example, the line starting with ‘#%cmd=’contains the specific command line parameters used when the report filewas generated.

Following the meta information is the summary metrics for each CEL file.The format is a column for each CEL file and a row for each metric witha tab separating each column. This file can be opened within Excel forinstance as a tab separated text file. The various quality controlmetrics reported are as follows with each metric potentially reportedmultiple times for different subsets of probesets. For example, the AFFXspike control probesets are reported separately (if included in theanalysis) and as part of the total.

The following terms are used herein and given the following meanings.Probe Count is the number of probes (features) analyzed. Probeset Countis the number of probesets analyzed. % Probes DABG Detected is thepercent of probes with a DABG probe level p-value less than or equal to0.01. This metric is only reported if the DABG quantification method isselected. % Probesets DABG Detected is the percent of probesets with aDABG probeset level p-value less than or equal to 0.01. This metric isonly reported if the DABG quantification method is selected. MeanProbeset PLIER Target Response: The mean PLIER signal estimate. MeanProbeset Abs PLIER Residuals is the mean absolute PLIER probe levelresiduals. Mean Probeset Abs PLIER RLE is the mean absolute PLIERrelative log expression (RLE). This metric is generated by taking thePLIER estimate for a given chip and calculating the difference in logbase 2 from the median value of that probeset over all the chips. Themean is then computed from the absolute RLE for all the probesets for agiven CEL file.

Positive vs Negative ROC AUC is the area under the curve (AUC) for areceiver operator curve (ROC) comparing the intron controls to the exoncontrols by applying a threshold to the PLIER signal estimate. The ROCcurve is generated by evaluating how well the PLIER signal estimateseparates the intron controls from the exon controls. The assumption(which is only valid in part) is that the intron controls are a measureof false positives and the exon controls are a measure of truepositives. An AUC of 1 reflects perfect separation whereas as an AUCvalue of 0.5 would reflect no separation. Note that the AUC of the ROCcurve is equivalent to a rank sum statistic used to test for differencesin the center of two distributions. These metrics in general (and theRLE, ROC AUC, and Residual metrics specifically) have proven most usefulin helping to identify outlier chip results relative to a given dataset.

Methods for deriving simple measures from array data, the feature or thefeature set level, to assess the quality of array data in terms ofcomparability are disclosed. It was argued that although comparingarrays in terms of feature intensity level is useful, decisionsregarding which arrays to exclude from analysis should be based onsignal level data—the data that will be used in downstream analysis.Other measures may also be used to assess quality.

Many of the quality assessment measures will be highly correlated. Forthe purpose of identifying gross outliers that should be excluded fromdownstream analysis any of the assessments discussed above will do agood job. Better methods—methods that are sensitive to departures inquality that have impact on particular downstream analyses—must by theirvery nature be developed on a case by case basis.

Conclusion

It is to be understood that the above description is intended to beillustrative and not do not by their details limit the scope of theclaims of the invention. Many variations of the invention will beapparent to those of skill in the art upon reviewing the abovedescription. The scope of the invention should be determined withreference to the appended claims, along with the full scope ofequivalents to which such claims are entitled. All cited references,including patent and non-patent literature, are incorporated herewith byreference in their entireties for all purposes.

1. A computer-implemented method for identifying exons that aredifferentially spliced between a first and second sample comprising: (a)obtaining an exon intensity measurement for a first exon in a first genein said first sample and in said second sample; (b) obtaining a genelevel measurement for said first gene in said first sample and in saidsecond sample, wherein the gene level measurement for said first sampleis a median normalized intensity value from one or more exons in saidfirst gene in said first sample and the gene level measurement for saidsecond sample is a median normalized intensity value from one or moreexons in said first gene in said second sample; (c) obtaining a firstnormalized intensity measurement for said first exon in said firstsample and a second normalized intensity measurement for said first exonin said second sample by dividing the intensity measurements obtained in(a) by the gene level measurements obtained in step (b) wherein the exonintensity in said first sample is divided by the gene level measurementin said first sample and the exon intensity in said second sample isdivided by the gene level measurement in said second sample; (d)calculating a splicing index measurement for said first exon, whereinthe splicing index is equal to the loge of the normalized intensity forsaid first exon in said first sample divided by the normalized intensityfor said first exon in said second sample; and (e) identifying saidfirst exon as being differentially spliced if the absolute value of thesplicing index obtained in (d) is greater than a threshold value.
 2. Themethod of claim 1 where the threshold value is 1-2.
 3. The method ofclaim 1 wherein the threshold value is
 2. 4. The method of claim 1wherein the normalized intensity measurements of step (c) are obtainedusing the following equation: $n_{i,j,k} = \frac{e_{i,j,k}}{g_{j,k}}$where e_(i,j,k) is the exon signal estimate of the i exon, j experiment,and k gene and g_(i,k) is the gene level signal estimate of the jexperiment and k gene.
 5. The method of claim 5 wherein a first value ofn_(i,j,k) is calculated for a first sample and a second value ofn_(i,j,k) is calculated for a second sample and the splicing index iscalculated by dividing the value of the n_(i,j,k) for said first sampleby the n_(i,j,k) for said second sample to obtain a third value andcalculating the log₂ of said third value.
 6. The method of claim 1wherein the exon is identified as being overrepresented in said firstsample relative to said second sample if the splicing index is greaterthan zero.
 7. The method of claim 1 wherein the exon is identified asbeing overrepresented in said first sample relative to said secondsample if the splicing index is greater than
 2. 8. The method of claim 1wherein the exon is identified as being underrepresented in said firstsample relative to said second sample if the splicing index is less than0.
 9. The method of claim 1 wherein the exon is identified as beingunderrepresented in said first sample relative to said second sample ifthe splicing index is less than −2.
 10. A computer-implemented methodfor identifying exons that are differentially spliced between a firstand second sample comprising: (a) obtaining an exon intensitymeasurement for a first exon in a first gene in said first sample and insaid second sample; (b) obtaining a gene level measurement for saidfirst gene in said first sample and in said second sample, wherein thegene level measurement for said first sample is a median normalizedintensity value from one or more exons in said first gene in said firstsample and the gene level measurement for said second sample is a mediannormalized intensity value from one or more exons in said first gene insaid second sample; (c) obtaining a first normalized intensitymeasurement for said first exon in said first sample and a secondnormalized intensity measurement for said first exon in said secondsample using the following equation:$n_{i,j,k} = \frac{e_{i,j,k}}{g_{j,k}}$ where e_(i,j,k) is the exonsignal estimate of the i exon, j experiment, and k gene and g_(i,k) isthe gene level signal estimate of the j experiment and k gene; (d)obtaining a splicing index measurement for said first exon, wherein afirst value of n_(i,j,k) is calculated for a first sample and a secondvalue of n_(i,j,k) is calculated for a second sample and the splicingindex is calculated by dividing the value of the n_(i,j,k) for saidfirst sample by the n_(i,j,k) for said second sample to obtain a thirdvalue and calculating the loge of said third value; and (e) identifyingsaid first exon as being differentially spliced if the absolute value ofthe splicing index obtained in (d) is greater than a threshold value.11. The method of claim 10 wherein the threshold value is between aboutzero and about
 2. 12. The method of claim 10 wherein the exon intensitymeasurements are obtained by hybridizing nucleic acid derived from thesample to an array of probes wherein there is a probeset comprising aplurality of probes for each exon to be interrogated and there are morethan 10,000 different probesets.
 13. A method for predicting if an exonin a gene is differentially spliced in a first and a second sample sets,comprising: obtaining a first exon level intensity measurement for saidexon in said first sample set and a second exon level intensitymeasurement for said exon in said second sample set; assume under a nullhypothesis of no alternative splicing for said exon that all exons ofsaid gene are expressed at levels that are proportional to each other insaid first and second sample sets; fit a model that predicts exonresponse under said null hypothesis; construct a statistic that measureshow deviant the observed data is from the model; use the statistic toconstruct a p-value; and analyze the p-value to predict if the exon isdifferentially spliced.
 14. The method of claim 13 wherein the model ise_(i,j,k)=α_(i,k)g_(j,k) where e_(i,j,k) is the signal of the i-th exonof the j-th sample of the k-th gene, g_(j,k) is the signal of k-th genein the j-th sample, and α_(i,k) is the ratio of exon i signal to itsgene signal.
 15. The method of claim 13 wherein the model ise_(i,j,k)=α_(j,k)p_(i,j,k)g_(j,k) where e_(i,j,k) is the signal of thei-th exon of the j-th sample of the k-th gene, g_(j,k) is the signal ofthe k-th gene in the j-th sample, α_(i,k) is the ratio of exon i signalto its gene signal in the sample where it is maximally expressed, and0≤p_(i,j,k)≤1 is a Splicing Index that estimates the proportionateexpression of this exon of this gene in tissue j.